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ABSTRACT 

The  Hydrographic  Sciences  Group  of  the  Naval  Postgraduate  School's 
Department  of  Oceanography  has  a  goal  to  establish  and  continue  a  research  base 
related  to  the  Global  Positioning  System  (GPS).  To  meet  this  goal  a  library  of 
programs  capable  of  reducing  GPS  satellite  data  to  position  and  time  solutions, 
independent  of  the  receiver's  navigation  solution,  is  being  acquired. 

This  thesis  describes  the  modification,  documentation,  and  establishment  of 
Texas  Instruments  4100  GEOSTAR  GPS  software  written  by  the  Naval  Surface 
Weapons  Center  for  their  CDC  Cyber  865  computer,  so  that  it  will  run  on  the  IBM 
3033  computer  at  NPS.  Position  and  time  solutions  for  either  a  static  or  dynamic 
receiver  can  be  achieved  using  this  software. 
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I.  INTRODUCTION 

The  Hydro  graphic  Sciences  Group  at  the  Naval  Postgraduate  School  (NPS), 
conducted  its  first  Sea  Floor  Benchmark  Experiment  in  May  1985.  The  experiment 
basically  was  to  deploy  a  benchmark  on  the  ocean  bottom  and  to  establish  its  position 
with  geodetic  accuracy.  This  was  to  be  accomplished  using  GPS  and  acoustic  range 
data  (Kumar,  et  al.,  1984). 

Beyond  the  goal  of  the  Sea  Floor  Benchmark  Experiment  in  May  1985  was  a 
higher  goal  of  establishing  and  continuing  a  GPS  research  base  at  NPS.  In  order  to 
reach  this  goal  it  has  become  evident  that  a  library  of  GPS  programs  should  be 
established  at  NPS  capable  of  reducing  GPS  data,  independent  of  outside  sources,  into 
a  format  compatible  for  analysis. 

The  Texas  Instruments  receiver  (TI-4100)  was  selected  as  one  of  the  GPS 
positioning  receivers  to  acquire  data  for  the  experiment.  One  version  of  the  TI-4100  is 
for  commercial  users  and  the  other  (TI-4100  GEOSTAR)  is  a  tri-agency  version  for  use 
by  the  Defense  Mapping  Agency  (DMA).  National  Oceanic  and  Atmospheric 
Administration  (NOAA).  and  the  United  States  Geological  Survey  (USGS).  TI-4100 
GEOSTAR  software  has  been  obtained  for  use  at  NPS  from  NOAA  and  the  Naval 
Surface  Weapons  Center  (NSWC)  in  Dahlgren,  Virginia,  for  compatibility 
modifications  with  NPS's  IBM  3033  mainframe  computer.  Of  this  software,  three 
programs  from  NSWC,  written  chiefly  by  S.  L.  Meyerhofffor  their  CDC  Cybcr-865 
computer,  reformat  GPS  data  and  then  compute  receiver  positions  and  time  by  either 
an  iterative  least  squares  approach  or  a  Kalman  filter. 

The  main  objectives  of  this  thesis  are  to  describe  the  effort  in  modifying  NSWC's 
software  codes  and  algorithms  for  compatibility,  and  establishing  them  on  the 
IBM-3033  mainframe  computer  at  NPS.  In  addition,  discussions  of  the  mathematical 
models,  results  of  test  runs  and  validation  of  the  modified  software  have  been  included. 
This  software  was  acquired  from  NSWC  for  exclusive  use  at  NPS.  Future  users  at 
NPS.  are  encouraged  to  test  and  verify  their  solutions  as  correct  before  accepting  them. 


II.  BACKGROUND 

A.       GLOBAL  POSITIONING  SYSTEM 

GPS  is  a  universal  satellite  positioning  system  that  is  independent  of  weather, 
time,  or  geographical  position.  It  provides  both  position  and  time  in  an  earth-fixed 
geocentric  cartesian  coordinate  system  as  well  as  velocity  data.  Since  it  is  a  space 
based  radio  navigation  system,  it  can  provide  accurate  information  to  an  unlimited 
number  of  users  globally  on  or  near  the  earth's  surface.  Three-dimensional  position 
data  to  within  a  16-m  spherical  error  probability  (SEP),  velocity  to  within  0.06-0.15  m/s 
and  GPS  system  time  to  within  25  ns  are  possible  instantaneously  to  properly  equipped 
users  anywhere  within  500  nmi  of  the  earth  (Milliken  and  Zoller,  1980). 

The  GPS  is  comprised  of  three  major  segments,  the  Space,  the  Control,  and  the 
User. 

1.    Space  Segment 

The  GPS  when  fully  operational  will  consist  of  18  active  satellites  (Space 
Vehicles  or  SVs)  in  six  planes  with  near  circular  orbits  at  altitudes  of  10,900  nmi  and 
periods  of  12  hours.  Each  of  the  six  orbital  planes  is  inclined  55  degrees  to  the  equator 
and  is  to  contain  three  SVs  equally  spaced.  The  configuration  of  the  SVs  in  the  orbital 
planes  is  such  that  a  minimum  of  four  SVs  will  be  available  to  the  user  at  all  times, 
anywhere  on  or  near  the  earth's  surface.  Each  individual  SV  transmits  its  broadcast 
ephemeris  and  clock  correction,  plus  an  almanac  which  contains  orbital  parameters  and 
clock  correction  estimates  for  all  other  SVs  in  the  system. 

Each  SV  transmits  two  radio  frequencies  simultaneously,  LI  and  L2 
respectively,  containing  navigation  information  in  the  navigation  data  message.  The 
LI  frequency  is  centered  on  1575.42  MHz  and  the  L2  frequency  is  centered  on  1227.60 
MHz.  These  signals  are  in  a  frequency  range  that  provides  good  all  weather  operation 
and  allows  the  user  to  determine  ionospheric  propagation  delay.  The  LI  signal  is 
modulated  by  a  10.23  MHz  clock  rate  Precision  (P)  code  and  a  1.023  MHz 
Course/Acquisition  (C/A)  code  (Spilkcr,  1980).  The  P-codc  is  a  random  sequence 
binary  code  that  provides  the  user  the  capability  for  high-precision  positioning  and  is 
highly  resistant  to  electronic  and  multipath  interference.  The  C/A-code  gives  all  users 
rapid  acquisition  capability  and  acts  as  an  aid  to  gain  access  to  the  P-code.    The 
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navigation  data  message  allows  the  user's  receiver  to  calculate  position,  velocity  and 
GPS  time  by  providing  broadcast  ephemeris,  clock  corrections  and  almanacs  for  all 
SVs.   The  L2  signal  is  modulated  by  the  P-code  only. 

2.  Control  Segment 

The  control  segment  consists  of  a  master  control  station,  a  number  of  monitor 
stations  and  ground  antennae  located  around  the  world.  The  USAF  master  control 
station  is  at  Colorado  Springs,  Co.,  and  monitor  stations  are  located  on  Kwajalcin 
Atoll,  Ascension  Island,  Hawaii  and  Diego  Garcia.  DMA  also  maintains  three  stations 
in  the  United  Kingdom,  Australia  and  Argentina,  and  stations  are  proposed  for 
Ecuador  and  Bahrain.  The  monitor  stations  acquire  ranging  data  for  each  SV  in  a 
passive  mode.  These  data  are  then  relayed  to  the  master  control  station  where  they 
are  processed  to  determine  the  SVs  position  and  signal  data  accuracy  (Wooden,  1985). 
The  master  control  station  takes  the  relayed  data,  computes  errors,  generates  a  new 
navigation  data  message  and  then  uploads  this  message  to  each  SV  via  the  ground 
antennae.  The  ground  antennae  transmit  as  well  as  receive  satellite  control 
information. 

3.  User  Segment 

This  segment  consists  of  stationary,  low,  medium  and  high  dynamic  receivers 
designed  with  specific  requirements  for  individual  users.  The  user  receivers  are 
designed  to  receive  and  process  SV  data  from  four  SVs  either  simultaneously  or 
sequentially.  The  TI-4100  receiver  determines  the  range  from  the  SV  by  measuring  the 
time  delay  of  the  SVs  specific  epoch  and  the  user  generated  codes.  Then  making 
separate  ranging  and  carrier  frequency  Doppler  shift  measurements  to  four  SVs,  the 
receiver  calculates  the  user's  position,  clock  bias  error,  clock  drift  and  velocity.  The 
position  coordinates  are  in  the  World  Geodetic  System  (WGS)  1972,  an  earth-fixed 
geocentric  cartesian  coordinate  system. 

B.       TEXAS  INSTRUMENTS  4100  NAVSTAR  RECEIVER  (TI-4100) 
1 .    Characteristics  and  Description 

The  TI-4100  receiver  is  portable  for  field  operations  (it  is  about  the  size  of  a 
small  suitcase)  but  can  be  mounted  in  a  standard  19-inch  rack  if  desired.  A  small 
lightweight  hand-held  Control  Display  Unit,  contains  a  keyboard  and  a  display 
window.  The  receiver  is  modular  in  design  and  allows  for  easy  removal  and 
replacement  of  circuit  boards. 
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The  TI-4100  is  a  single  channel  dual-frequency,  multiplexing  NAVSTAR 
receiver  and  navigation  processor  that  tracks  up  to  four  SVs  (Texas  Instruments, 
1983).  The  dual  frequency  receiver  allows  for  tracking  LI  and  L2  frequencies  for  each 
of  the  four  SVs.  It  also  has  an  optional  dual  cassette  recorder  for  recording  GPS/SV 
data.  If  GPS/SV  data  is  to  be  acquired  in  a  degraded  mode  (i.e.  less  than  four  SVs) 
there  is  the  capability  of  adding  an  external  atomic  clock  time  standard  input  to  the 
TI-4100  receiver.  The  TI-4100  receiver,  Control  Display  Unit,  antenna  and  cassette 
recorder  are  shown  in  Figure  2.1. 

2.  Receiver  Measurements  and  Postion  Calculations 

The  basic  measurements  made  by  the  TI-4100  receiver  are  the  arrival  times  of 
the  SV-generated  signals  compared  to  the  receiver  clock,  which  contains  a  stable 
oscillator  and  a  counter  (Texas  Instruments,  1983).  Four  pseudorange  measurements 
plus  the  navigation  data  message  provide  enough  information  to  calculate  a 
three-dimensional  user  position  solution  plus  the  user  clock  correction  to  true  GPS 
time.  The  solutions  are  obtained  using  four  observation  equations  and  four  unknowns. 
Doppler  shifts  due  to  the  relative  velocity  between  the  SV  and  the  user  are  measured. 
Measurements  are  also  made  of  the  phase  differences  for  the  LI  and  L2  signals 
between  the  SV  and  the  user. 

3.  Position  Accuracy 

Error  statistics  for  the  TI-4100  receiver  are  given  in  Table  I.  It  should  be 
noted  for  this  table  that  Operating  Modes  using  time  aiding,  have  a  positive  growth 
rate  error  of  2.0  to  10.0  m/h  (1-sigma)  for  frequency  standards  that  allow  navigation 
calculations  to  estimate  precise  frequency  (Texas  Instruments,  1983).  The  positioning 
accuracies  in  Table  I  were  based  on  the  TI-4100  receiver  observing  pseudoranges, 
Doppler  shift  data,  phase  differences  for  LI  and  L2,  then  computing  position,  time  and 
velocity  solutions.  The  pseudoranges  from  the  TI-4100  receivers  are  the  major 
observables  considered  in  this  thesis. 

C.       PSEUDORANGE  DATA  AND  THEIR  CHARACTERISTICS 

The  basic  operation  of  the  GPS  depends  on  the  user  determining  pseudorange 
and  range  rate  to  a  number  of  GPS  SVs  which  have  precisely  known  ephemerides 
(Texas  Instruments,  1983).  Pseudorange  is  the  apparent  range,  which  includes  the 
clock  errors  between  the  S  V  and  the  user's  receiver.  The  transit  time  of  the  navigation 
signal  between  each  SV  and  the  user  is  measured  by  the  receiver  and  then  scaled  by  the 
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Figure  2.1     TI-4100  NAVSTAll  System 
(From  Texas  Instruments,  I  1-4100  Owner's  Manual,  19S3). 
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TABLE  I 

TI-4100  POSITION  ACCURACY 

(From  Texas  Instruments,  TI-4100  Owner's 

Manual,  1983) 

1-Sigma 
Position  Error 

1-Sigma 
Velocity 

Operating  Mode 

Error 
Stati  sties 

P 

(m) 

C/A 

Error 
(m/s) 

Four  satellites 

Stationary 

0. 5  g  acceleration 

2  g  acceleration 

4  q  acceleration 

Differential 

Stationary  differential 

3-D 
3-D 
3-D 
3-D 
3-D 
3-D 

13 
14 
14 
30 
4 
3 

46 
47 
47 
55 
12 
8 

N/A 

0.15 

2.00 

5.00 

0.15 

N/A 

Three  satellites 

Degraded  (altitude  aiding) 

Horizontal 

16 

24 

0.15 

Differential  degraded 
(altitude  aiding) 

Horizontal 

5 

6 

0.  15 

Degraded  (external  cesium 
frequency  standard  time 

3-D 

16 

47 

0.15 

aiding) 

Differential  Degraded 

(external  cesium  frequency 
standard  time  aiding) 

3-D 

5 

12 

0.15 

Two  satellites 

Degraded  (altitude  and 
time  aiding) 

Horizontal 

16 

24 

0.15 

Differential  degraded 

(altitude  and  time  aiding) 

Horizontal 

5 

6 

0.15 

speed  of  light  to  compute  a  raw  pseudorange.  Each  raw  pseudorange  must  be 
corrected  for  tropospheric  and  ionospheric  propagation  delays,  SV  clock  oiTset,  user 
clock  bias,  relativity  corrections,  and  earth  rotation  corrections.  The  pseudorange 
calculation  is  given  by  Equation  (2.1). 

R*  -  R.+  cAtAi  +  c(Atu-Atsi)  (2.1) 

where  i  =   1  to  4  (pertaining  to  four  observed  SV's) 

14 


R-*    =  pseudorange  to  the  SV 

R:     =  true  range 

c      =  the  speed  of  light 

At<jj  =  SV^  clock  offset  from  GPS  system  time 

Atu    =  user  clock  ofiset  from  GPS  system  time 

^t-Ai  =  propagation  delays  and  other  errors 

The  tropospheric  delay  is  estimated  using  simple  refraction  models  with  standard 
atmospheric  values.  Ionospheric  delay  is  calculated  by  taking  advantage  of  the  dual 
frequency  measurements  of  the  LI  and  L2  signals.  The  SV  clock  error  (difference 
between  SV  clock  and  GPS  system  time)  is  modeled  by  the  Control  Segment  and  these 
correction  coefficients  are  transmitted  in  the  navigation  data  message  from  the  SV  to 
the  user.  User  clock  bias  is  due  to  the  fact  that  the  receiver  clock  is  not  perfectly 
synchronized  to  the  SV  clock.  Relativity  corrections  are  due  to  the  eccentricity  of  the 
SV's  orbit.  Earth  rotation  corrections  are  due  to  the  motion  of  the  earth  during  the 
signal's  propagation  from  SV  to  the  user.  Once  four  raw  pseudoranges  have  been 
corrected  for  systematic  errors,  they  are  used  to  calculate  the  receiver's 
three-dimensional  position  and  time. 

The  precision  of  position  calculated  by  pseudoranges  is  dependent  on  the 
geometry  of  the  SV's  with  respect  to  the  user.  This  geometry  is  continually  changing, 
even  if  the  user's  receiver  is  not  in  motion  relative  to  the  earth.  If  the  intersection  of 
the  pseudoranges  is  orthogonal,  errors  in  the  pseudorange  measurement  are  minimized 
whereas  if  the  angle  of  intersection  is  small,  errors  are  large  because  of  the  geometric 
relationship. 

Geometric  Dilution  of  Precision  (GDOP)  is  a  composite  "indicator"  expressing 
the  influence  of  the  SV's  geometry  on  the  accuracy  of  the  user's  clock  ofiset  and 
position  determinations  (Milliken  and  Zoller,  1980).  Other  terms  to  describe  dilution 
of  precision  in  position  and  time  are:  Position  Dilution  of  Precision  (PDOP)  which 
yields  dilution  of  precision  in  three-dimensions;  Horizontal  Dilution  of  Precision 
(HDOP)  which  gives  dilution  of  precision  in  the  two  horizontal  dimensions;  Vertical 
Dilution  of  Precision  (VDOP)  which  yields  the  dilution  of  precision  in  the  vertical 
dimension;  and  finally  Time  Dilution  of  Precision  (TDOP)  is  the  dilution  of  precision  in 
the  user's  time  bias  (Milliken  and  Zoller,  1980).  GDOP,  which  is  dimcnsionless,  is 
calculated  by  liquation  (2.2),  (Milliken  and  Zoller,  1980), 
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GDOP  =  [(PDOP)2  +  (TDOP)2]1/2  (2.2) 

The  best  geometry  for  a  four  SV  configuration  is  the  combination  which  results  in  the 
lowest  GDOP. 

Once  pseudorange,  Doppler  shift  and  phase  data  have  been  acquired,  and  the 
receiver's  field  solutions  calculated,  post-processing  of  these  data  may  yield  improved 
position  and  time  solutions.  To  make  TI-4100  data  accessible  to  other  users  with 
different  computer  processing  systems,  it  has  become  necessary  to  format  SV  data  in  a 
manner  that  can  be  exchanged  among  independent  users. 

D.       STANDARDIZED      EXCHANGE       FORMAT      FOR      NAVSTAR      GPS 
GEODETIC  DATA 

This  Standardized  Exchange  Format  (SEF)  was  proposed  by  V.  Dan  Scott  and  J. 
Gary  Peters  (1983),  and  was  adopted  to  make  possible  the  universal  exchange  of 
NAVSTAR  GPS  data.  To  do  this,  serious  compromises  were  made  in  regard  to 
optimum  magnetic  tape  utilization,  computer  execution  efficiency,  and  computer 
program  size  (Scott  and  Peters,  1983).  These  compromises  were  required  to  allow  for 
maximum  tape  compatibilty  between  different  computer  systems. 

1.    The  Basic  Format 

The  building  block  for  each  standardized  data  exchange  tape  is  the  record, 
which  is  80  characters.  Under  this  format,  an  unlabeled  9-track,  ASCII-coded  tape 
with  an  80-character  blocking  factor  is  used.  The  header  (or  the  first)  file  on  a  tape 
consists  of  comment  records  and  is  not  considered  part  of  the  data  set  but  is  used  to 
help  the  user  to  identify  the  contents  of  the  tape.  When  the  last  data  file  is  written  on 
the  tape  a  special  end  of  information  file  is  written  (Scott  and  Peters,  1983). 

On  the  tape,  records  fall  into  two  categories:  Control  Information  or  Data 
Information  records.  Control  Information  records  are  used  to  identify  which 
FORTRAN  Format  statement  to  use  when  reading  data  records.  Data  Information 
records  consist  of  one  of  eight  data  items  that  are  uniquely  identified  by  a  four  digit 
data  item  number.   The  eight  data  items  under  the  SEF  are  listed  below: 

(1)  Constant  equipment  and  software  specific  data 

(2)  Campaign  and  station  specific  data 

(3)  Calibration  data 

(4)  Time  tag  update  data 

(5)  Space  Vehicle  specific  data 
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(6)  Weather  data 

(7)  Real  time  solution  data 

(8)  Measurement  data 

The  goal  of  this  format  is  to  be  independent  of  hardware  and  software  design.    In 
addition  to  data  formats,  the  reference  Coordinate  System  for  the  data  must  be  defined. 
2.    Coordinate  System 

The  WGS  72  an  earth-fixed  geocentric  cartesian  coordinate  system  is  defined 
as  follows  (Seppelin,  1974): 

X-Axis  =  Intersection  of  the  WGS  72  reference  meridian  plane 
and  the  plane  of  the  Equator. 

Y-Axis  =  Measured  in  the  plane  of  the  Equator  90  degrees 
east  of  the  X-Axis. 

Z-Axis  =  Parallel  to  the  direction  of  the  conventional 

international  origin  (CIO)  for  polar  motion. 

The  reference  ellipsoid  for  WGS  72  is  defined  by  the  parameters  semi-major  axis,  a  = 
6,378,135  m,  and  flattening,  f  =  1  298.26  (Seppelin,  1974). 
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III.  DATA  PROCESSING 

The  pre-processing  of  GPS  data  is  carried  out  to  obtain  a  format  suitable  for 
further  manipulation.  This  is  necessary  when  geodetic  positioning  of  the  receiver 
(either  static  or  dynamic)  is  desired,  through  post-processing  the  SV  data.  Figure  3.1 
portrays  the  flowchart  for  processing  the  TI-4100  data. 

A.      PRE-PROCESSING 

1.  Manipulation  of  Cassette  Recorded  Data 

The  first  step  in  processing  GPS  data  recorded  on  cassettes,  is  to  dump  the 
raw  data  from  the  cassettes  to  either  a  9-track  magnetic  tape  or  disk.  The  Applied 
Research  Laboratories  (ARL)  at  Austin,  Texas  have  developed  software  to  accomplish 
this  fundamental  task.  The  software  at  ARL  can  dump  and  process  the  raw  data 
recorded  on  cassettes,  and  format  the  GPS  data  onto  a  9-track  magnetic  tape  in  the 
SEF  (Chapter  H.D.).  This  tape  can  then  be  transferred  to  disk  for  convenience  in 
further  processing  with  a  mainframe  computer. 

2.  Alignment  of  Common  Time-Tagged  Data 

The  alignment  of  common  time-tagged  data  items  is  accomplished  through  use 
of  the  CON9TR  program. 

a.    Program  Description 

Program  CON9TR,  the  First  of  three  programs  acquired  from  NSWC 
(Chapter  I.),  reads  the  SEF  data  and  formats  it  into  another  NSWC  specified  File.  This 
file  then  becomes  the  input  for  the  two  remaining  programs,  SOLTOM  and  KALMN2 
and  has  all  common  time-tagged  data  on  the  same  record.  The  only  calculations 
performed  at  this  stage  are  variances  on  the  pseudorange  data  and  the  delta  range 
(Dopplcr)  data  which  are  measures  of  the  signal  to  noise  ratio. 

The  program  also  looks  at  the  "quality  vectors"  (Texas  Instruments, 
1982)  of  the  data  on  the  9-track  SEP"  tape  and  returns  an  integer  value  representative 
of  the  data  quality  at  that  time  line.    Four  integer  values  are  possible: 

0  =  good  data 

1  =  good  data  but  this  is  the  first  Doppler  count 

14  =  bad  data  at  that  point  but  there  was  not  a  break 

in  the  Doppler  data 

15  =  bad  data 

18 


I 1-4100 

CASSETTES 
GPS  DATA 


9 -TRACK    STANOARO 
EXCHANGE    FORMAT  TAPE 
FOR   GPS   DATA 


USER    INPUT 
FILE 


PROGRAM  CON9TR 


MESSAGE 
FILE 


USER    INPUT 
FILE 
(SOLUTION  TYPE) 


STATIC   I    DYNAMIC 


PROGRAM 
SOUOM 


RESIOUA1 
FILE 


PROGRAM 
KAl  MN2 


GO  OP 

FItE 

(OPTIONAL) 


RECEIVER  POSITION 

AY   SOLUTION 

TYPE   ANO    TIME 


USER    INPUT 
FILE 
(SOLUTION   TYPE) 


POSITION 
HISTORY 


RECEIVER   POSIT  1014 
BY   SOLUTION 
TYPE    AND    TIME 


Figure  3.1     Flowchart  For  Processing  TI-4I00  GEOSTAR  GPS  data. 
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The  CON9TR  program  requires  a  user  input  file  (Figure  3.1)  in  addition 
to  the  GPS  data  in  the  Standardized  Exchange  Format  during  program  execution. 
CON9TR  consists  of  the  main  program  and  22  subroutines.  Appendix  A  gives  a  brief 
description  of  each  subroutine's  function  and  nesting. 

b.  Input/Output 

The  user  input  file  (Figure  3.1)  for  the  CON9TR  program  consists  of  the 
following: 

(1)  The  number  of  trackers,  usually  four  but  may  be  less  than  four 

(2)  GPS  time  of  the  first  data  point  to  be  put  on  the  new  file 

(3)  Time  to  end  data  on  the  new  file,  if  the  end  of  the  file 

is  found  before  the  time  to  end  data,  then  data  will  end 
at  the  end  of  file 

(4)  Range  variance  bias  factor  [km  ] 

(5)  Doppler  variance  bias  factor  [km  ] 

(6)  X-coordinate  of  the  initial  receiver  position  that 

will  be  put  on  the  new  files 

(7)  Y-coordinate  of  the  initial  receiver  position  that 

will  be  put  on  the  new  files 

(8)  Z-coordinate  of  the  initial  receiver  postion  that 

will  br  put  on  the  new  files 

Appendix  B  gives  the  user  input  guide  to  CON9TR  program  variables  and  format. 

The  first  output  file  from  CON9TR  (Figure  3.1)  writes  all  common  time 
satellite  tracking  data  into  the  same  record  block,  and  Appendix  C  gives  the  output 
format.  The  second  output  file  writes  out  messages  to  aid  the  user  in  data  quality 
analysis,  flagging  one  of  four  quality  vector  conditions  that  exist  in  the  data. 

c.  Compution  of  Variances 

The  variances  for  the  pseudorange  and  Doppler  data  in  the  CON9TR 
program  are  computed  by  Equations  (3.1)  and  (3.2),  respectively  (Texas  Instruments, 
1982). 

a2code  =  KBL  *  N)/(10(Ci/10))]  *  (3.1) 

[0.5  +  BpD/(10(Ci/10))]  *  Kj  +  K2 

where  the  asterick  ,*,  is  used  to  signify  multiplication 
B^  =  code  loop  noise  bandwidth  Hz 

20 


BpQ  =  code  loop  predetection  bandwidth  Hz 

C-  =  L-  carrier  power  to  noise  spectral  density  ratio 

dB-Hz 
N  =  tracking  mode,  1,  2,  3,  4  based  on  SV  mode 
Kj  =  (chips)"  to  [m~]  conversion: 

(29.305)2  for  P-code 

(293.05)2  for  C  A-code 
K->  =  variance  bias  factor  [nV] 


(T2pLL  =  [(BL  *  N)  (10(Ci  10))1  *  {1  +  [0.5  *  (3.2) 

BpD(10,Ci  10))]}  *K3  -  K4 
where, 

PLL  =  carrier  phase-locked  loop 
B^  =  carrier  loop  bandwidth  Hz 
^PD  =  carrier  predetection  bandwidth  Hz 
K;  =  [rad~J  to  [m~]  conversion 
=  (c  27rf)2 
c  =  speed  of  light 
f  =  L  band  carrier  frequency 
Lj  =  1.57524  x  109  Hz 
L^  =  1.2276  x  109  Hz 
Kj  =  variance  bias  factor  [m^] 

d.    Program  Modifications  at  NPS 

The  original  code  from  NSWC  was  written  in  FORTRAN  IV.  This 
code  was  modified  at  NPS  to  be  compatible  with  the  FORTRAN  77  compilers 
(FORTVS  and  WF77)  on  the  IBM-3033  mainframe  computer.  Changes  to  the  original 
code  included  the  following. 

•  Changes  in  the  formatted  write  statements 

•  Changes  in  variable  alignments  in  the  common  statements 

(i.e.  making  sure  double  precision  variables  are  aligned 
ahead  of  single  precision  variables) 

•  Changes  in  dimension  size  of  the  data  arrays 

•  Changes  in  defining  variables 
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•  Changes  in  the  placement  of  data  intialization 

statements 

•  Addition  of  the  Block  Data  subroutine  to  contain  the 

moved  data  intialization  statements  and  accompanying 
common  statements. 

Executive  files,  written  at  NPS  to  handle  the  input  and  output  files  for 
the  CON9TR  program,  included  the  auto-double  precision  commands  for  all  real 
variables  using  the  FORTVS  compiler.  This  modification  automatically  doubles  the 
declared  precision  in  the  code  and  was  required  due  to  the  use  of  60  bit  words  by  the 
CDC  Cyber-865  computer  and  the  32  bit  words  by  the  IBM-3033  mainframe  computer 
at  NPS. 

B.       DETERMINATION  OF  PSEUDORANGE  POSITION  SOLUTION 

User  receiver  position  solutions  whether  fixed  (static)  or  moving  (dynamic)  may 

be  calculated  by  post-processing  the  TI-4100  GEOSTAR  GPS  data. 
1.    Static  Positioning 

The  second  program,  SOLTOM  from  NSWC  (Chapter  I.)  computes  a  position 

solution  for  a  static  receiver  by  using  an  iterative  batch  least  squares  calculation. 

a.  Program  Description 

Program  SOLTOM  computes  receiver  position  and  time,  solving  for 
initial  corrections  to  estimates  of  parameters  X,  Y,  Z,  time  biases  and  scaling  factor  for 
a  tropospheric  correction.  Types  of  solutions  which  can  be  selected  through  the  user 
input  file  are  described  in  Chapter  III. B. I.e..  The  program  consists  of  the  main 
program  and  52  subroutines.  Appendix  D  gives  a  brief  description  of  each  subroutine's 
function  and  nesting. 

b.  Math  Models 

The  math  models  for  the  batch  least  squares  solution  in  program 
SOLTOM  solve  for  corrections  to  the  initial  estimates  of  the  parameters  desired.  The 
observation  equation  is  given  by  Equation  (3.3): 

R0,  =  [(Xsv-Xp)2  +  (YSV-  Yp)2  +  (Zgy-Zpft^  +  CB  +  CD  +  CA  +  TR       (3.3) 
where, 


R°-  =  i  n  range  observed  from  SV  of  receiver 
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Xsv  =  earth-fixed  geocentric  X-coordinate  of  SV 
Ysv  =  earth-fixed  geocentric  Y-coordinate  of  SV 
Z$v  =  earth-fixed  geocentric  Z-coordinate  of  SV 

X     =  earth-fixed  geocentric  X-coordinate  of  receiver 

Y     =  earth-fixed  geocentric  Y-coordinate  of  receiver 

Z     =  earth-fixed  geocentric  Z-coordinate  of  receiver 

CB  =  Clock  bias  term,  SV  or  receiver  depending  on  user 
input  selection 

CD  =  Clock  drift  term,  SV  or  receiver  depending  on  user 
input  selection 

CA  =  Clock  aging  term  (frequency  drift),  SV  or  receiver 
depending  on  user  input  selection 

TR  =  Tropospheric  refraction  correction. 

The  least  squares  solution  of  Equation  (3.3)  in  matrix  form,  can  be  written  as: 

BP  +  E  =  0  (3.4) 

where, 

B  =  Matrix  containing  the  partial  derivatives  of  the  observation 
equation  with  respect  to  the  parameters  being  solved  for 

P  =  Vector  containing  the  differences  of  the  parameter  adjusted 
(Pa)  and  the  parameter  estimate  (Pe),  (Pa  -  Pe) 

E  =  Vector  containing  the  partial  derivatives  of  the  observation 
equations  multiplied  by  the  weight  matrix,  multiplied  by  the 
data  residuals 

The  matrix  B  is  obtained  from  Equation  (3.5), 

B  =  ATwA  (3.5) 
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where, 

A  =  Partial  derivatives  of  data  equation  with  respect  to  a 
parameter  j 

A     =  Transpose  of  the  A-matrix 

w  =  Weight  matrix  with  c\  on  the  diagonal,  zeros 

off  diagonal,  where  k  =   1  to  number  of  data  points,  this 
is  the  variance  of  noise  level  on  a  "k"  range  measurement 

or  writing  Equation  3.5  in  summation  form, 

Bi,i  =  Ink  =  i  ^ROk  /  ^  *  5R\  /  apj  *  V  0-«> 

where, 

n  =  number  of  data  points  used  in  the  least  squares  solution 

R^o  =  ktn  data  equation 

P^  =  i     parameter  being  solved  for 

P:  =  jin  parameter  being  solved  for 

th 

Gy2  =  k     variance  value 
The  A-matrix  or  design  matrix  is  defined  by  Equation  (3.7): 

Ai,j  =  dR°i  I  5Pj      '  (3-7) 

where, 

R°i  =  ith  data  point  (3.3) 

P:  =  j"-"  parameter  being  solved  for 

The  inverse  of  B  is  the  variance-covariance  matrix,  COV,  with  the  variances  (p-  of  the 
parameter  corrections  AP  that  are  being  solved  for,  on  the  diagonal  given  by  Equation 
(3.8), 
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B"1  =  COV 


where.  COY  is  a  6  x  6  matrix  if  six  parameters  are  being  solved,  or  7  x  7  if  the  optional 
1  parameter  •  CA  is  to  be  determined.  The  off  diagonal  elements  are  the 
cross-correlations  of  the  paired  parameters  for  that  time  line.  The  vector  E  in 
summation  form  is  shown  in  Equation  (3.9), 

Ei-Znk=l^Rk0^pi»*ffk2*RreS  <3-9> 

where. 

Rjl  a  =  ktb  data  equation 

o  =  number  of  data  poir.'.s  used  in  Least  squares 
calculation 

rh 
P:  =  i     parameter  being  solved  for 

2  i 

cr-.~  =  i"'  -.ve-.z.-.-.ec  varia:  .  - 

RTes  =  range  residuals 

or  writing  Equation  (3.9)  in  matrix  form, 

E  =  Rrs  -.vA  5.10) 

""■here. 

A?  =  B-1E  3.11> 

ir.h 

AP  =  corrections  to  the  initial  estimates  of  unknown  parameters  being  solved  for   X. 
Y,  Z.  dock  biases,  and  tropospheric  refraction  correction). 

The  following  Equations  (3.12)  to  (3.35)  describe  the  various  corrections,  Time  Epoch, 
clock,  Ionospheric  refraction,  Tropospheric  refraction,  earth  rotation,  and  rclativistic 
effects  which  are  required  to  correct  pscudorange  and  delta  range  (Doppler)  data 
-Meyerhoff,  19_;  5 


Time  Epoch  Correction 


Tsv=tt-Toc  (3.12) 

where, 
Tgv  =  time  from  epoch 

tt  =  the  GPS  time-tag  for  range  and  delta  range  (Doppler) 
data 

TQC  =  time  of  reference  for  clock  corrections 

To  correct  the  T  v  for  the  end  of  the  week  cross  over,  correct  Tsv  as  follows: 

If  Tsv  <  302400  [s]  then,  Tsv  =  Tsv  +  604800  [s] 
If  Tsv  >  302400  [s]  then,  Tsv  =  Tsv  -  604800  [s] 

Clock  Correction 


Cd(tt)  =  [A0  +  (A!  *  Tsv)  +  (A2  *  Tsv  *  Tsv)]  *  c  (3.13) 


ACd(tt)  =  [A0  +  (Aj  *  {(2  *  Tsv)  -  td})]  *  td  *  c  (3.14) 

where, 

C  i(tt)  =  clock  correction  for  pseudorange  data 
at  GPS  time-tag  in  [km] 

ACc|(tt)  =  clock  correction  for  delta  range 

(Doppler)  data  at  GPS  time-tag  in  [km] 

AQ  =  clock  correction  bias  term  [s] 

A  j  =  clock  correction  drift  term  [s/s] 


A2  =  clock  correction  aging  term 
[s/(s)2] 

T$v  =  time  from  epoch 
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c  =  speed  of  light  299792.458  [km/s] 

td  =  the  Doppler  count  interval  for  delta  range 
(Doppler)  data 

Ionospheric  Refraction  Correction 

Ci0(tt)  =  CPj(tt)  -  CP2(tt)/(Q1/Q2)2  -  1.0  (3.15) 


DCi0(tt)  =  [(DOPT  ^tO/Q,)  -  (DOPT  ?(tt)/Q?)]  *  c  (3.16) 

[(Qi/Q2)2  -  i.O]  vos 

where, 
CjQ(tt)  =  pseudorange  ionospheric  refraction 
delay  in  [km] 

tt  =  GPS  time-tag  for  pseudorange  and  delta  range 
(Doppler) 

CPj  =  LI  pseudorange  in  [km] 

CP2  =  L2  pseudorange  in  [km] 

DOPj^j  =  LI  Doppler  count 

DOPj^2  =  L2  Doppler  count 

DCj0  =  delta  range  (Doppler)  frequency  shift  due 
to  ionosphere 

Qj  =  LI  frequency  multiplier  (154.0) 

Q2  =  L2  frequency  multiplier  (120.0) 

VQS  =  The  nominal  satellite  oscillator  frequency,  10.23  MHz 

In  Equation  (3.15)  the  correction,  C-  (tt),  is  used  to  correct  the  pseudorange  that  is 
time-tagged,  tt.  In  Equation  (3.16)  the  correction,  DCj0(tt),  is  applied  to  the  delta 
range  (Doppler)  data  to  correct  it  for  time-tag,  tt. 
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Tropospheric  Refraction  Correction 

The  Tropospheric  refraction  computation  can  be  divided  into  two  parts, 
weather  dependent  and  receiver  geometry  calculations  (Meyerhoff,  1985). 
Weather  Dependent  Correction 

First,  convert  temperature  to  degrees  Kelvin  and  convert  humidity  to  a 
fraction: 

Tk  -  Tc  +  273.0  (3.17) 


RH  =  Hd/ 100.0  (3.18) 

where, 
T    =  air  temperature  in  degrees  centigrade 

H^  =  the  relative  humidity  in  percent. 

Second,   compute   the   surface   partial   pressure,    EQ,   from   temperature   and   relative 
humidity. 

E0  =  RH  *  35.65  *  10  exp[7.617-(2285.0/Tk)]  (3.19) 

Third,  the  components  of  the  zenith  value  of  tropospheric  delay,  Z^      and  Zwet  are 
calculated  from,  EQ,  and  atmospheric  pressure. 

Zdry  =  2.276  *  Pb  *  0.01  (3.20) 


Zwet  =  (470  *  E0exp(1.23)/Tk)  +  (1.705  *  106  (3.21) 


Ah*E0exp(1.46/Tk3) 


where, 
exp  =  exponent 


Pb  =  the  atmospheric  pressure  in  [mb] 
A^  =  temperature  lapse  rate  set  at  0.006  °C/m 
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exp  =   exponent  It  should  be  noted  that  Zd      and  Zwet  need  only  be  recomputed, 
when  the  weather  data  changes. 

Receiver  Geometry  Dependent  Correction 

First,  to  find  the  elevation  angle,  EV(t),  of  the  satellite  at  time,  t,  take 
the  inverse  sine  of  the  dot  product  of  the  unit  vector  from  the  receiver  to  the  satellite 
and  the  unit  position  vector  of  the  receiver. 

EV(t)  =  sin-1    RS  ~R  (3.22) 

MGrs  '   MGr 
where, 
RS  =  vector  from  the  receiver  to  the  satellite 

MGrs  =  vector's  magnitude 


R  =  position  vector  of  the  receiver 
MGr  =  magnitude  of  R 

MGrs  =  [(X(t)  -  Xr)2  +  (Y(t)  -  Yr)2  4-  (Z(t)  -  Zj.)2]1/2  (3.23) 

MGr  =  [Xr2  +  Yr2  +  Zj.2]1'2  (3.24) 

where, 

X  ,  Y  ,  Z    =  the  estimate  of  the  receiver's 

position  in  an  earth-fixed 

geocentric  coordinate  system 

rewriting  Equation  (3.23),  EV(t)  can  be  written  as, 

SEV(t)  =  [{(X(t)  -  Xr)  *  Xr)  (3.25) 

+  {(Y(t)  -  Yr)  *  Yr)} 
+  {(Z(t)  -  Zr)  *  Zr}]/(MGrs  *  MGr) 
where, 
X(t),  Y(t),  Z(t)  =  the  position  coordinates  in  a  carth- 
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fixed  geocentric  coordinate  system 
or 

EV(t)  =  sin-1  [SEV(t)]  (3.26) 

Second,  the  multiplication  factors   for  dry  and  wet  as  a  function  of  elevation  are 
calculated  from  the  elevation  angle,  EV,at  GPS  time-tag  (tt). 

Fdry(t)  =  1  /   [sin  EV(t)J  +  {0.00143/[(tan  EV(t))+  0.0445]}  (3.27) 


Fwet(t)  =  1  /   [sin  EV(t)]  +  {0.00035/[(tan  EV(t))+  0.17]}  (3.28) 

Third,  Fdp,(t  -  td)  and  Fwet(t  -  td)  are  calculated  in  a  similar  manner  as  Equations 
(3.27  and  3.28). 

Total  Tropospheric  Refraction  Correction 
TR(tt)  =  [Zdry  *  Fdry(t)]  +  [Zwet  *  Fwet(t)]  (3.29) 

DTR(tt)  =  {Zdry  *  [Fdry(t)  -  Fdry(t  -  td)]}  +  (Zwet  *  [Fwet(t)  -  (3.30) 

W1  -  td)» 

where, 

TR(tt)  and  DTR(tt)  are  range  corrections  are  in  [km]. 

Earth  Rotation  Correction 

An  earth  rotation  correction  is  required  due  to  the  rotation  of  the  earth 
during  the  signal  propagation. 

Ccr  =  [(Y(t)  *  Xr)  -  (X(t)  *  Yr)]  *  RRE/c  (3.31) 
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DCer  =  [({Y(t)  -  Y(t  -  td)}  *  Xr)  -  ({X(t)-X(t  -  td)}  *  Yr)]  *  RRE/c  (3.32) 

where, 

X(t),  Y(t),  Z(t)  =  position  coordinates  in  an  earth- 
fixed  geocentric  cartesian 
coordinate  system  at  pseudorange 
transmit  time  "t"  used  with  data 
time  tagged  "tt" 

Xf,  Yr      =  estimate  of  the  X  and  Y  coordinate 
in  an  earth-fixed  geocentric 
cartesian  coordinate  system 
RRE  =  Rotational  rate  of  the  Earth  corrected  for  motion 
of  the  equinox 
(0.00007292115855  [rads/s] 

td  =  Doppler  count  interval  for  delta  range 
(Doppler)  data. 

Relativistic  Correction 

The  correction  for  the  relativistic  effects  is  due  to  the  eccentricity  of  the 
satellite's  orbit. 

Rel(tt)  =  (2  *  [GM]1'2/^  *  e  *  A  *  sin  [Ek(t)]  (3.33) 

at  the  beginning  of  the  Doppler  count, 

DRel(tt)  =  (2  *  [GM]'/2/c)  *  e  *  A  *  sin  [Ek(t  -  td)]  (3.34) 

which  yields, 

DRC(tt)  =  RC(t)  -  RC(t  -td)  (3.35) 

where, 

tt  =  GPS  time-tag  for  both  pseudorange  and  delta  range 
(Doppler)  data 
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td  =  Dopplcr  count  interval  for  delta  range  (Doppler) 
data 

t  =  GPS  time-tag  at  transmission  time  for  data  time- 
tagged  tt 

e  =  eccentricity  of  the  satellite  orbit 

A  =  the  square  root  of  the  sejni-major  axis  of  the 
satellite's  orbit 

Ek(t)  =  eccentric  anomaly  for  the  satellite  at  pseudorange 
transmission  time  t 

Ek(t  -  td)=  eccentric  anomaly  for  the  satellite  at  the 
beginning  of  the  Doppler  count  interval 

GM  =  WGS  72  value  of  the  Earth's  gravitational  constant  including 
atmosphere  398600.8  km3/s2 

RC  =  computed  range 

these  are  the  defined  constants  and  variables. 

c.  Input/Output 

There  are  two  input  files  for  the  Soltom  program  (Figure  3.1).  The  first 
file  is  the  user  input  file,  where  the  user  sets  which  types  of  solutions  will  be  calculated 
and  which  data  solutions  will  be  printed  out.  The  user  input  guide  for  this  file  is 
described  in  Appendix  E.  The  second  input  file  for  SOLTOM  is  an  output  file  from 
the  CON9TR  program  containing  records  that  are  aligned  with  common  time-tagged 
GPS  data  items. 

There  are  also  two  output  files  for  SOLTOM  (Figure  3.1).  The  first  file 
is  the  updated  receiver  position,  time  and  other  print  options  such  as  solution  sigmas 
and  clock  bias  terms,  (as  determined  by  the  user  input  file).  The  second  output  file  is  a 
data  corrections  residual  file,  which  can  be  used  in  plotting  routines  and  for  input  to 
other  software  that  handle  multiple  station  positioning  (McyerhofT,  1986). 

d.  Program  Modifications  at  NPS 

The  source  code  from  NSWC  for  the  SOLTOM  program  was  written  in  FORTRAN 
IV  which  has  been  modified  to  run  at  NPS  on  the  IBM-3033  mainframe  computer. 
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Similar  changes  to  the  source  code  in  the  CON9TR  program  were  also  made  to  the 
SOLTOM  program  including  the  executive  files  (Chapter  1 1 1. A. 2.).  The  original  code 
was  written  to  be  able  to  handle  the  precise  ephemeris  data  as  well  as  the  onboard 
broadcast  ephemeris  data  in  the  navigation  data  message  from  the  SV.  The  test  data 
for  SOLTOM  was  generated  at  NSWC  using  the  Broadcast  ephemeris  only,  and  as 
such,  the  modified  program  version  at  NPS  has  the  Precise  Ephemeris  subroutines  and 
corresponding,  call,  statements  commented  out. 

e.    Solution  Techniques  Through  Least  Squares 

There  are  ten  types  of  solutions  possible  with  the  SOLTOM  program 
which  are  dependent  on  the  user  input  file  (Meyerhoff,  1985).  Appendix  E  gives  the 
details  for  the  type  of  solution  desired.  Of  these  ten  types,  nine  are  batch  least  squares 
solutions  and  one  solution  type  five  is  a  Newton-Raphson  solution  (a  geometrical 
solution,  not  a  least  squares  fit).  The  type-five  solution  can  be  used  for  either  static  or 
dynamic  positioning  performing  an  independent  solution  at  each  time  line.  The  nine 
batch  least  squares  solution  types  fall  into  three  data  classes,  Range  (pseudorange), 
Range  Difference  (delta  range  Doppler),  and  Bias  Range. 

The  reader  is  WARNED,  that  the  Bias  Range  (solution  type  ten)  has  still 
not  been  tested  by  NSWC,  and  thus  it  should  not  be  used  at  NPS.  Solution  types  one, 
four,  six,  and  nine  are  Range  (pseudorange)  solutions.  Solution  types  two  and  seven 
are  Range  Difference  (Doppler)  solutions.  Finally,  solution  types  three  and  eight  are 
combination  solutions  of  the  Range  (pseudorange)  and  Range  Difference  (Doppler) 
data. 

The  SOLTOM  program  has  the  capability  to  run  more  than  one 
iteration  of  batch  least  squares,  if  needed  (e.g.  for  bad  initial  receiver  position  or  time 
bias  problems).  The  user  has  control  via  the  user  input  file  to  set  how  often  the  least 
squares  solution  will  be  calculated  and  printed  out.  When  a  solution  is  to  be 
computed,  the  program  takes  the  B-matrix  containing  the  partial  derivatives  of  data 
equations  with  respect  to  the  parameters  being  solved  for.  and  multiplies  by  the 
variances  of  the  range  measurements  (Equations  3.5  and  3.6),  inverts  the  B-matrix, 
which  is  the  variancc-covariance  matrix  with  variances  of  the  corrections  AP  on  the 
diagonal,  and  computes  the  corrections  to  the  initial  estimates  of  the  parameters  being 
sought,  AP  (Equation  3.11).  Then  Equation  (3.36)  gives  the  updated  station  position 
coordinates,  clock  bias,  drift,  aging,  and  Tropospheric  refraction  correction: 


X  =  X  +  APX  (3.36) 

33 


Y  =  Y  +  APV 
> 


Z  =  Z  +  APZ 


Ao  =  Ao  +  APAo 
A,t  =  Ajt  +  APAjt 


A2t2  =  A2t2  4-  APA 


TR  =  TR  +  APTR 

where, 

X,  Y,  Z  =  updated  earth-fixed  geocentric  cartesian 
coordinates  through  batch  least  squares 

AQ  =  clock  correction  bias  term  in  [s] 

Aj  -  clock  correction  drift  term  in  [s/s] 

A2  =  clock  correction  aging  term  in 
[s/(s)2] 
(optional  7th  parameter  to  solve  for) 

TR  =  tropospheric  refraction  correction 

The  new  solution  takes  all  previous  solutions  and  averages  up  to  the  current  time  line, 
before  producing  the  iterated  solution. 
2.    Dynamic  Positioning 

The  third  program,  KALMN2  from  NSVVC  computes  a  position  solution  for  a 
moving  (dynamic)  receiver,  using  the  UDU"  factorization  of  an  eight-state  extended 
Kalman  filter  (Chapter  I.). 

a.    Program  Description 

Program  KALMN2  uses  the  first  four-states  for  receiver  position  (X,  Y, 
Z)  and  receiver  clock  bias,  and  the  other  four-states  for  receiver  velocity  along  three 
axes  (X,  Y,  Z)  and  receiver  clock  drift.  By  using  the  Doppler  phase  (change  in  range) 
measurements  to  smooth  the  pseudoranges,  signal  multipath  and  measurement  noise 
effects  can  be  suppressed  significantly  to  produce  a  viable  solution  (Mcycrhoff  and 
Evans,  1986). 
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Note:  The  four-states  of  the  velocity  and  clock  drift  are  still  under 
testing  at  NSWC. 

The  KALMN2  program  consists  of  the  main  program  and  43 
subroutines.  Appendix  F,  gives  a  brief  description  of  each  subroutine's  function  and 
nesting. 

b.    Math  Models 

The  math  model  for  the  UDU'  factorization  of  the  inverse  of  the 
B-matrix  (which  contains  the  partial  derivatives  of  the  observation  equations  with 
respect  to  the  parameters  being  solved  for)  used  in  program  KALMN2  is  described  in 
Bierman  (1977)  and  Meyerhoff(1986). 

The  TI-4100  receiver  provides  two  pseudorange  and  two  Doppler 
measurements  (via  LI  and  L2  frequencies),  for  as  many  as  four  satellites 
simultaneously.  The  following  Equations  (3.37  to  3.42)  describe  the  various  corrections 
and  conversions;  receiver  clock  bias,  continuous  count  (Doppler)  to  interval  (Doppler) 
data,  and  Ionospheric  refraction,  which  are  required  for  the  Doppler  and  the 
pseudorange  data  (for  the  Ionospheric  correction)  (MeyerhofT,  1985): 

Receiver  Clock  Bias  Correction 

tt  =  ttr  -  RCB  (3.37) 

where, 
ttr  =  receiver  clocks  GPS  time-tag  for  range  and  delta 
range  (Doppler)  data 

tt  =  GPS  time-tag  corrected  for  receiver  clock  bias 

RCB  =  estimate  of  the  receiver  clock  bias  which  is 

initially  set  to  zero  and  an  improved  estimate 
of  "RCB"  is  made  after  each  time  line 

The  time  delay  data  are  multiplied  by  the  speed  of  light  "c"  and 
converted  to  a  pseudorange  in  [km]. 

Continuous  Count  to  Interval  {Doppler)  Data  Conversion 

DOPuCtt^  =  DC^tt^  -  DC^tt  -  td)i  -  (td  *  OFFSET1)  (3.38) 
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DOPL2(tt)i  =  DC2(tt)i  -  DC2(tt  -  td)i  -  (td  *  0FFSET2)  (3.39) 

where, 
td  =  Doppler  count  interval  for  delta  range  (Doppler)  data 

DCjCtt  -  td)  ,  DC2(tt  -  td)=  LI  and  L2  continuous 

Doppler  count  data  at  the  delta  range  interval 

DCi(tt)  ,  DC2(tt)  =  LI  and  L2  continuous  count  Doppler 

data  at  the  end  of  the  delta  range  interval 
OFFSET j  ,  OFFSET2  =  Doppler  offsets  used  by  the  receiver 

for  LI  and  L2,  -6000  and  +7600  respectively 

DOPL1(tt)  ,  DOPL2(tt)  =  LI  and  L2  interval  Doppler  data 
for  time-tag  "tt" 

i  =  satellite  from  which  data  originates 

The  delta  range  (Doppler)  data  is  obtained  by  converting  the  LI  interval  Doppler  data 
for  each  satellite. 

RD(tt)i  =  (-DOPLH(tt)/V0S  *  Qj)  *  c  (3.40) 

where, 
VQS  =  nominal  satellite  oscillator  frequency,  10.23  MHz 

Q j  =  LI  frequency  multiplier  (154) 

RD(tt)  =  delta  range  (Doppler)  for  one  satellite  in  [km] 
Ionospheric  Refraction  Correction 

RCI(tt)i  =  CPIOty  +  ClOftOj  (3.41) 


RDCI(tt)i  =  RD(tt)i  -  DCIO(tt)i  (3.42) 

where, 
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CPI  =  uncorrected  pseudorange  measurement  in  [km] 

RD(tt)  =  uncorrected  delta  (Doppler)  range  measurement 

ClO(tt)  =  ionospheric  refraction  correction  in  [km] 

for  pseudorange  data  (refer  to  Equation  3.15) 

DCIO(tt)  =  ionospheric  refraction  corrections  in  [km] 

for  delta  (Doppler)  .range  data  (refer  to  Equation 
3.16). 

Smoothing  For  Noise 

The  pseudorange  data  are  contaminated  due  to  signal  multipath  and 
other  noise  sources  (Meyerhoff  and  Evans,  1986).  To  suppress  the  noise  due  to 
multipath  in  the  pseudorange  data,  they  are  smoothed  by  the  delta  range  (Doppler) 
data.   The  next  set  of  Equations  (3.43,  3.44,  and  3.45)  describes  this  smoothing: 

BN(tt)  =  R(tt)  -  DR(tt.ttfst)  (3.43) 


Bavg(")  =  IBj(u)/N  (3-44) 


Rsm(tt)  =  Bavg(tt)  +  DR(tt,ttfst)  (3.45) 

where, 

tt  =  GPS  time-tag  for  pseudorange  and  delta  range  (Doppler)  data 

ttfst  =  GPS  time-tag  for  the  first  good  pseudorange  and 
Doppler  count  of  this  data  span 

R(tt)  =  Ionospheric  corrected  pseudorange  in  [km] 

DR(tt, ttfst)  =  Ionospheric  corrected  delta  range  (Doppler)  data 

N  =  number  of  data  points  in  the  span 

B^-  =  N     estimate  of  pseudorange  at  time  "ttfst" 
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B„v„  =  average  value  of  the  first  pseudorange 

6 

R     (tt)  =  smooth  pseudorange,  if  the  absolute  value  of 

the  difference  between  the  raw  pseudorange  and 
smoothed  pseudorange  exceeds  20  [m],  the 
data  point  is  considered  bad. 

Both  the  pseudorange  and  delta  range  (Doppler)  data  must  be  corrected  for  Time 
Epoch,  clock  errors,  earth  rotation,  relativistic  effects,  and  Tropospheric  refraction. 
These  formulae  are  the  same  ones  used  for  the  iterative  least  squares  solution  in  the 
SOLTOM  program  (Equations  3.12,  3.13,  3.14,  3.31,  3.32,  3.33,  3.34,  3.35,  3.29,  3.30). 
The  corrected  pseudorange  is  given  by, 

RCCtty  =  R^tty  +  Cd(tt)i  -  RelCtt);  -  [TR^  *  (1  +  CR)]  +  Cer  (3.46) 

RDC(tt)  is  calculated  in  a  similar  manner  as  (3.46)  for  the  delta  range  (Doppler)  data. 

The  following  Equations  (3.47  to  3.51)  describe  computed  ranges  and  residuals  for 
pseudorange  and  delta  range  (Doppler)  data  (Meyerhoff,  1985): 

Computed  Ranges 

CRCtt);  =  [(X(t)i  -  Xr)2  +  (Y(t)rYr)2  4-  (Z(t\  -  Zr)2]1/2  (3.47) 

for  the  range  at  the  beginning  of  the  Doppler  interval, 

CR(tt  -  td)j  =  [{X(t2)j  -  [Xr  -  (XVr  *  td)]}2  +  {Y(t2)j  -  [Yr  -  (YVr  *  td)]}2    (3.48) 

+  {Z(t2)j  -  [Zr  -  (ZVr  *  td)]}  2]1/2 

where, 

X(t),  Y(t),  Z(t)  =  position  of  ith  satellite  in  an  earth- 
fixed  geocentric  cartesian  coordinate 
system  at  transit  time  t 

Xr,  Yr,  Zr  =  estimate  of  the  receiver's  position  in  an 

earth-fixed  geocentric  cartesian  coordinate 
system 
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X(t2),  Y(t2),  Z(t2)  =  position  of  ith  satellite  in  an 

earth-fixed  geocentric 
cartesian  coordinate  system  at 
time  t2 

XVr,  YVr,  ZVr  =  estimate  of  receiver's  velocity 

tt  =  GPS  time-tag  corrected  for  receiver  clock,  bias 

td  =  Doppler  count  interval  for  delta  range  (Doppler) 
data 

Therefore  the  computed  range  difference  for  the  individual  satellites  is, 

CRDftty  =  CR(tt)j  -  CR(tt  -  td)i  (3.49) 

Residuals 

The  pseudorange  residuals  are  given  by, 

RomcftOj  =  RC(tt)j  -  CR(tt)i  -  RCB  (3.50) 

and  for  the  delta  range  (Doppler)  data, 

RDomc(tt)i  =  RDC(tt)j  -  CRD(tt)i  -  RCD  (3.51) 

where, 
RC(tt)-  =  corrected  pseudorange  (Equation  3.46) 

RDC(tt)-  =  corrected  delta  range  (Doppler) 

CR(tt):  =  computed  pseudorange 

DCR(tt)j  =  computed  delta  range  (Doppler) 

RCB  =  estimate  of  the  receiver  clock  bias 

RCD  =  estimate  of  frequency  bias  of  the  receiver  oscillator 
which  is  also  receiver  clock  drift. 
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tt  =  GPS  time-tag  corrected  for  receiver  clock  bias 

td  =  Doppler  count  interval  for  delta  range  (Doppler) 
data 

i  =  i      satellite  being  tracked 

A  state  vector  SP  is  corrected  and  improved  over  time,  is  composed  of 
the  receiver  position  in  an  earth-fixed  geocentric  cartesian  coordinate  system  (Xr,  Yr, 
Zr)  and  receiver  clock  bias  (RCB).  The  last  set  of  Equations  (3.52  to  3.58)  describes 
the  observation  equation  and  development  of  the  state  vector  for  the  pseudorange  data 
(Meyerhoff,  1985): 

Observation  Equation 

D(u){  =  CRCtty  +  RCB  (3.52) 

or  rewriting  it  as, 

D(tt)i  =  [(X(t\  -  Xr)2  +  (Y(t)i  -  Yr)2  +  (Z(t)x  -  Zr)2]1/2  +  RCB  (3.53) 

The  H-matrix  is  defined  as  containing  the  partial  derivatives  of  the 
observation  equation  with  respect  to  the  state  vector: 

H^j  =  aD(tt)/5Xr  =  -(X(t)i  -  Xr)/CR(tt)i  (3.54) 

U2-  -  0D(tt)/dYr  =  -CY(t)j  -  Yr)/CR(tt)i 

H3'j  =  dD(u)/dZr  =  -(Z(t)j  -  ZrVCRCtOj 

H4i  =  aD(tt)/5RCB  =  1 

The  PA-matrix  is  defined  as  a  4x4  variance-covariance  matrix  for 
position  and  receiver  clock  bias,  with  oil  diagonal  elements  zeroed  out  and  the  diagonal 
elements  set  to  the  initial  variance  of  Xr,  Yr,  Zr,  and  RCB.  The  Kalman  vector  K  is 
given  by  the  matrix  equation, 

K(tt)j  =  PA(tt-)  *  HTi  *  [  Hj  *  PA(tt-)  *  HTj  +  Rj  J'1  (3.55) 

where, 
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PA(tt-)  =  variance-covariance  matrix  at  time  tt  before  it  is 
updated  for  this  data  point 

R     =  variance  on  the  pseudorange  observation 

i     =  i      satellite  being  tracked 

Once  the  PA-matrix  is  updated  it  is  defined  by, 

PA(tt+)  =  PA(tt-)  -  K(tt)j  *  H{  *  PA(tt-)  (3.56) 


State  Vector 

The  final  two  equations  define  the  updated  state  vector  SP(tt  +  ): 

SP(tt  +  )  =  SP(tt-)  +  K(tt)i  *  Romc(tt)i  (3.57) 


or 


Xr(  +  )  =  Xr(-)  +  KjCtt^  *  Romc(tt)i  (3.58) 

Yr(  +  )  =  Yr(-)  +  K2(tt)i  *  Romc(tt)1 
Zr(  +  )  =  Zr(-)  +  K3(tt)i  *  Romc(tt)i 
RCB(  +  )  =  RCB(-)  +  K4(tt)j  *  Romc(tt)i 

where  (  +  )  means  the  updated  state  vector  and  (-)  means  before  the  updated  state 
vector. 

c.    Input  Output 

There  are  two  input  files  for  the  KALMN2  program  (Figure  3.1).  The 
first  file  is  an  output  file  from  the  CON9TR  program  which  contains  records  that  are 
aligned  with  common  time-tagged  GPS  data  items.  The  second  file  is  a  user  input  file 
which  sets  the  types  of  solutions  and/or  constraints  on  that  solution,  along  with  which 
data  solutions  will  be  printed  out.   The  user  input  guide  for  this  file  is  in  (Appendix  G). 

There  are  three  possible  output  files  for  KALMN2  (Figure  3.1).  The 
first  file  is  the  updated  receiver  position  and  time  and  any  other  print  options  that  arc 
set  by  the  user  input  file.  The  second  file  is  a  position  history  file  that  is  used  in 
printing  out  the  updated  receiver  position  and  time  file.  The  third  file  is  an  optional 
GDOP  file,  that  can  be  set  in  the  user  input  file. 
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d.  Program  Modifications  at  NPS 

The  source  code  from  NSWC  for  the  KALMN2  program  was  written  in 
FORTRAN  IV  which  has  been  modified  to  run  at  NPS  on  the  IBM-3033  mainframe 
computer.  Similar  changes  to  the  source  code  for  the  CON9TR  and  SOLTOM 
program  were  also  made  to  the  KALMN2  program  including  the  executive  files 
(Chapter  III.A.2.).  The  original  code  was  written  to  handle  both  the  precise  as  well  as 
the  onboard  broadcast  ephemeris  data  in  the  navigation  data  message  from  the  SV. 
The  test  data  for  KALMN2  was  generated  at  NSWC  using  the  broadcast  ephemeris 
only,  and  as  such,  the  modified  program  version  at  NPS  has  the  precise  ephemeris 
subroutines  and  corresponding,  call,  statements  commented  out. 

e.  Solution  Techniques  Through  Kalman  Filter 

Using  a  UDU'  factorization  of  an  eight-state  extended  Kalman  filter 
(with  smoothed  pseudorange  measurements)  in  program  KALMN2,  there  are  nine 
types  of  solutions  that  may  be  calculated  depending  on  the  user  input  file  (Appendix 
G). 

Note  that,  as  was  stated  previously  in  Chapter  III. B. 2. a.,  the  velocity 
solutions  from  the  four-velocity  states  are  still  being  tested  at  NSWC. 

The  initial  receiver  position  (X,  Y,  Z)  and  clock  biases  are  updated  using 
the  corrected  state  vector  from  the  Kalman  filter  (Equations  3.57  and  3.58).  If  for  a 
given  time  line  the  absolute  value  for  the  ionospheric  correction  ClO(tt)  exceeds  20  m, 
that  data  is  considered  bad  and  flagged  as  such,  with  a  quality  vector  set  to  15.  Height 
constrained  (especially  in  marine  areas)  allows  the  user's  receiver  to  track  four  satellites 
with  poor  GDOP  and  still  be  able  to  produce  a  good  solution. 
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IV.  TEST  RUN  RESULTS  AND  VALIDATION 

The  testing  of  the  three  programs  CON9TR,  SOLTOM,  KALMN2  was 
accomplished  by  using  two  independent  data  sets  in  the  SEF  on  9-track  magnetic  tape 
as  received  from  NSWC.  Included  also  were  hardcopy  printouts,  listing  the  solutions 
from  those  data  sets  using  the  CON9TR,  SOLTOM,  and  KALMN2  programs  on 
NSWC's  CDC  Cyber-865  computer.  A  direct  comparison  of  the  output  listing  files 
containing  the  receiver  position  solutions  and  times  for  the  SOLTOM  and  KALMN2 
programs  from  both  the  NSWC's  CDC  Cyber-865  and  the  NPS  IBM-3033  computers 
could  then  be  made. 

After  the  original  source  codes  and  some  algorithms  were  modified  for 
compatibility  on  the  NPS  IBM-3033,  the  programs  were  test  run,  using  the  FORTVS 
compiler  with  an  auto-double  precision  option  discussed  previously  (Chapter  1 1  LA. 2.). 

A.  STATIC  POSITION  SOLUTION  RESULTS 

Table  II  shows  the  differences  between  results  computed  at  NPS  and  NSWC 
using  the  NSWC  test  data  set  (from  Henderson  Point,  85006)  for  the  SOLTOM 
program.  Running  the  CON9TR  and  SOLTOM  programs,  position  solution  results 
(X,  Y,  Z)  were  duplicated  to  the  millimeter  level  for  solution  type  seven;  to  the 
centimeter  level  for  solution  types  one  and  two;  to  the  decimeter  level  for  solution 
types  four  and  six;  and  to  the  meter  level  for  solution  type  nine  (Appendix  E). 

The  above  differences  in  the  SOLTOM  solutions  are  considered  insignificant  due 
to  the  limitations  of  the  observed  data  (Meyerhoff,  1986). 

B.  DYNAMIC  POSITION  SOLUTION  RESULTS 

Table  III  shows  the  differences  between  results  (beginning  and  ending  solutions 
using  the  test  data  set)  computed  at  NPS  and  NSWC  using  the  NSWC  test  data  (from 
Henderson  Point,  95005)  for  the  KALMN2  program.  Running  the  CON9TR  and 
KALMN2  programs,  position  solutions  and  clock  biases  were  duplicated  to  the 
millimeter  and  10      second  level  respectively. 

It  should  be  noted  that  not  all  options  (Chapter  III.B.2  and  Appendix  G)  within 
the  programs  have  been  checked  and  tested.  At  this  point  it  can  be  said  that  these 
programs  have  been  established  and  validated  on  the  IBM-3033  at  NPS  to  the  degree 
of  being  able  duplicate  NSWC's  results  if  the  same  user  options  are  followed. 
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TABLE  II 

CON9TR  AND  SOLTOM  SOLUTION  RESULTS 

D  iff  ere 

nee  Between 

NPS  and  NSWC  Solutions 

X.Y.Z 

Solution       Position 
Type       (kilometers) 

Clock  Bias 
(seconds) 

CI 

ock  Drift 
(seconds) 

Of 

Number 
Data  Points 

ONE 
Range  Data 
Receiver  Clock 

10"5 

lO"9 

10"11 

1133 

TWO 
Doppler  Data 
Receiver  Clock 

10'5 

N/A 

10"11 

1125 

FOUR 
Common   Range 
Receiver  Clock 

10-4 

N/A 

N/A 

832 

SIX 

Range  Data 
SV  Clock 

lO"4 

lO"7 

10-" 

1133 

SEVEN 
Doppler  Data 
SV  Clock 

10"7 

N/A 

10"11 

1125 

NINE 

Common   Ranqe 
SV  Clock 

10"3 

10"6 

10"12 

832 
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TABLE  III 

CON9TR  AND  KALMN2  SOLUTION  RESULTS 

Difference  Between  NPS  and  NSWC  Solutions 

GPS  Time-Tag 
(seconds) 

(X.Y.Z)                Clock  Bias 
(kilometers)        (seconds) 

Number  of 
Doppler  Edited 
Points 

First  Solution 

14129.00 

10'5                        10"9 

N/A 

Last  Solution 

17499.00 

10"5                        10'6 

9 

45 


V.  SUMMARY 

I  low  to  process  the  GPS  satellite  data  from  raw  cassettes  off  the  TI-4100 
GEOSTAR  receiver  to  a  reduced  position  and  time  solution  has  been  discussed. 

The  capability  of  processing  raw  data  from  cassettes  to  SEF  at  NPS  is  a  task  for 
the  future. 

NPS  at  this  time  has  the  capability,  starting  with  a  SEF  9-track  magnetic  tape,  to 
use  the  CON9TR  program  (where  the  common  time-tagged  data  are  aligned),  to  select 
either  the  SOLTOM  (an  iterative  least  squares  for  positioning  a  static  receiver)  or  the 
KALMN2  program  (which  uses  a  Kalman  filter  for  positioning  a  dynamic  receiver)  and 
to  compute  position  and  time  solutions.  The  calculations  for  program  CON9TR,  and 
the  math  models  for  programs  SOLTOM  and  KALMN2  were  documented.  The 
changes  required  to  modify  the  original  codes  and  algorithms  for  these  programs 
written  for  NSWC's  CDC  Cyber-865  computer  to  make  them  compatible  with  the 
IBM-3033  computer  at  NPS  have  been  described.  Testing  of  the  programs  shows  a 
successful  conversion  (with  the  restrictions  discussed  in  Chapter  IV.)  of  the  programs 
to  the  IBM-3033  mainframe  computer.  This  software  forms  a  basis  for  a  GPS 
programs  library  that  will  aid  research  using  the  GPS  at  NPS. 
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VI.  RECOMMENDATIONS 

Since  the  conversion  of  the  CON9TR,  SOLTOM,  and  KALMN2  programs  has 
been  successfully  tested  using  the  two  independent  data  sets  from  NSWC,  it  is 
recommended  to  explore  other  options  within  the  programs  set  by  the  user  input  files 
for  the  SOLTOM  and  KALMN2  programs.  Options  that  should  be  explored  for  the 
SOLTOM  program  include  the  use  of  the  precise  ephemcris  and  clock  files,  the  use  of 
more  than  one  iteration  for  the  solutions,  and  the  calculation  of  Newton-Raphson 
range  solutions  followed  by  a  comparison  with  the  iterative  batch  least  squares 
solutions  for  a  given  time  line.  The  KALMN2  program  should  be  explored  by  using 
the  precise  ephemeris  and  clock  files,  turning  the  height  constraint  on  and  off  in  the 
user  input  file  and  comparing  the  solution  results.  The  refraction  correction  should  be 
turned  on  and  off  to  see  how  the  solutions  are  affected  (this  is  another  option  in  the 
user  input  file). 

It  is  further  recommended  to  explore  the  possiblities  of  rewriting  or  adding  codes 
and  algorithms  to  improve  the  solution  corrections  for  updating  the  four-state  vectors 
of  the  velocity  in  the  extended  eight-state  Kalman  filter  making  them  reliable  for  use  in 
velocity  solutions. 

In  the  future  TI-4100  software  for  carrier  phase  data,  differential  positioning  etc., 
should  be  acquired  and  modified  to  be  compatible  with  the  IBM-3033  computer  at 
NTS. 
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APPENDIX  A 
CON9TR  SUBROUTINE  DESCRIPTIONS  AND  NESTING 


INPUT  -  This  routine  reads  the  user  inputs. 

GETDAT  -  This  routine  reads  the  next  set  of  satellite 

tracking  data  from  the  9-track  Standard  tape 
and  also  updates  weather,  Navigation  data, 
almanac  and  calibration. 

RDDTL   -  This  routine  initializes  the  quality  vector 

flags  and  controls  branching  of  the  program 
depending  on  a  particular  error  Hag. 

ER9TR  -  This  routine  controls  branching  of  the  program 

for  several  error  and  terminating  conditions. 

UDNTDT  -  This  routine  processes  non-tracking  data  so  that 

it  can  be  used,  after  it  determines  which  type 
of  non-tracking  data  was  read  from  the  9-track 
Standard  tape. 

PRTSL   -  This  routine  prints  tape  solutions. 

UDWTDT  -  This  routine  gets  the  weather  data  from  the  9-track 

Standard  tape  and  stores  it  in  common  that  will  be 
used  to  calculate  the  tropospheric  correction. 

UDSATD  -  This  routine  reads  the  9:track  Standard  tape  and 

defines  new  satellite  orbit  clock  parameters. 

FDTRN   -  This  routine  will  give  the  tracker  number  for  the 

satellite  that  the  new  data  is  for. 

DEDPOF  -  This  routine  stores  the  values  for  the  LI  and  L2 

Doppler  offsets  so  that  they  may  be  used  in  the 
program. 

PRTPDT  -  This  routine  prints  any  unused  9-track  Standard 

tape  data. 

UDTRD   -  This  routine  decides  if  the  tracking  data  is 

continuous  count  Doppler  range  or  data  quality 
measurement  blocks  and  then  calls  the  correct 
subroutine  to  handle  the  correct  data  type. 

DEQUDT  -  This  routine  defines  data  qualitv  and  then  stores 

the  data  quality  information.    The  satellite 
tracker  mode  and  the  quality  vectors  are  used  to 
determine  if  the  data  are  good.    The  remaining  data 
values  in  this  block  are  used  to  calculate  errors 
on  the  tracking  data. 

DERGDT  -  This  routine  stores  the  range  data  and  the  signal 

to  noise  data  in  correct  data  array  that  will  Be 
processed. 

DEDPDT  -  This  routine  converts  the  accumulated  Doppler 

phase  data  into  range  difference  data.    The 
Doppler  phase  data  is  corrected  for  offset  and 
then  it  is  differenced  with  preceding  data  if  the 
preceding  data  was  good. 
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GNINP  - 

GNID    - 

IN9TR   - 

BLDATA - 
SDAT    - 

FSIG     - 

BLOCK  DATA 


This  routine  reads  the  logical  unit  number  of  the 
input  file  and  returns  the'information  read  one 
item  at  a  time. 

This  routine  controls  how  data  items  are  read  by 
reading  a  lour  digit  data  type  number. 

This  routine  initializes  the  read  package  of  the 
9-track  Standard  tape  and  sets  up  a  data  block 
for  the  data  items  read  of!  of  the  tape. 

This  routine  contains  the  data  that  has  been  read 
off  of  the  9-track  Standard  tape. 

This  routine  takes  the  site  name  and  the  data 
interval  from  the  9-track  Standard  tape  read 
common  and  passes  these  values. 

This  routine  determines  the  difference  between 
the  raw  range  and  the  range  difference  data  for 
the  satellite^being  tracked." 

This  routine  was  added  at  NPS  to  make  this 
program  compatible  with  the  FORTRAN  77 
compilers  on  the  IBM-3033  mainframe  computer. 
It  consists  of  common  statements  and  data 
statements  that  the  rigorous  compilers 
at  NPS  would  not  perinit  where  thev  were  placed 
in  the  original  source  code  of  CON9TR. 


MAIN  PROGRAM  CON9TR 
INPUT 
GETDAT 

RDDTL 

ER9TR 
UDNTDT 

PRTSL 

LDWTDT 

LDSATD 

FDTRN 
DFDPOF 
PRTPDT 
LDTRD 

DEQUDT 

FDTRN 
DERGDT 

FDTRN 
DEDPDT 

FDTRN 
GNINP 

GNID 
IN9TR 

BLDATA 
GNINP 
RDDTL 
SDAT 
FSIG 
BLOCK  DATA 
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APPENDIX  B 
INPUT  GUIDE  FOR  PROGRAM  CON9TR 

(  by  S.  L.  MeyerhofT,  NSWC) 

This  is  the  user  input  file  for  the  program  with  variable  names  and  descriptions.  The 
input  is  format  free. 

number  of  trackers  (4) 

GPS  time  of  first  data  point  to  put  on 
new  file 

time  to  end  data  on  new  file,  if  EOF 
found  before  THS  then  data  will  end  at 
EOF 

range  variance  bias  factor  (kilometer) 

doppler  variance  bias  factor  (kilometer) 


X  coordinate  of  initial  receiver  position 
that  will  be  put  on  the  new  files 

Y  coordinate  of  initial  receiver  position 
that  will  be  put  on  the  new  files 

Z  coordinate  of  initial  receiver  position 
that  will  be  put  on  the  new  files 


line  1 
1: 

NSAT    = 

2: 

TSS      = 

3: 

TES      = 

4: 

K2 

5: 

K4 

line  2 

1: 

X 

2: 

Y 

3: 

Z 
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APPENDIX  C 
TEXAS  INSTRUMENTS  GPS  SATELLITE  RECEIVER  FILE  FORMAT 

(  by  S.  L.  Meyerhoff,  1984,  NSWC  ) 
Data  from  the  DMA  version  of  the  TI-4100  GPS  satellite  receiver  can  be  put 
into  the  following  format.    In  this  format  all  common  time  satellite  tracking  data  will 
be  in  the  same  record  block.    The  format  for  this  file  is  as  follows. 


var 


Record  One 


word  1    I  TYPE 


Record  Two 


type 


intg 


meaning 


This  is  the  data  type  for 
the  next  record 


If  I  TYPE  equals  one  then  record  two  will  contain  initial  information. 


word  1  RLIOFF 

word  2  RL20FF 

word  3  PB 

word  4  TP 

word  5  MD 

word  6  X 

word  7   Y 

word  8   Z 

word  9   NTR 

word  10  DTREC 

word  11  SITE(l) 
word  12  SITE(2) 


real 
real 
real 
real 
real 
real 

real 

real 

intg 

real 

char 
char 


frequency  bias  offset  for  LI 

frequency  bias  offset  for  L2 

pressure  in  (millibars) 

temperature  (degrees  Celsius) 

relative  humidity  (percent) 

X  coordinate  of  initial 

receiver  position  (kilometers) 

Y  coordinate  of  initial 

receiver  position  (kilometers) 

Z  coordinate  of  initial 

receiver  position  (kilometers) 

maximum  number  of  trackers  used 

will  be  4  or  larger 

user  given  data  interval 

(seconds) 

site  identification  for  data 

site  identification  for  data 
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If  ITYPE  equals  two 
data  message) 

word  1   J 
word  2    ISAT(J) 
word  3   AS(1) 
word  4  AS(2) 

word  5  AS(3) 

word  6  CICS 

word  7   CISS 

word  8   CRCS 

word  9   CRSS 

word  10  CUCS 

word  11  CUSS 

word  12  DNS 

word  13  ES 
word  14  I  DOTS 

word  15  I  OS 
word  16  MOS 
word  17  0MEDS 

word  18  OMEGS 
word  19  SQAS 


(continued) 
then  record  two  will  contain  satellite  ephemeris  data  (navigation 


intg 

tracker  number 

intg 

satellite  prn  number 

real 

satellite  clock  bias  in  seconds 

real 

satellite  clock  drift 

(seconds/second) 

real 

satellite  clock  aging 

(seconds/second*second) 

real 

cross  track  cosine  amplitude 

(radians) 

real 

cross  track  sine  amplitude 

(radians) 

real 

radial  cosine  amplitude 

(kilometers) 

real 

radial  sine  amplitude 

(kilometers) 

real 

in  track  cosine  amplitude 

(radians) 

real 

in  track  sine  amplitude 

(radians) 

real 

corrections  to  mean  motion 

(radians) 

real 

eccentricity 

real 

inclination  rate 

(radians/second) 

real 

inclination  (radians) 

real 

mean  anomaly  (radians) 

real 

right  ascension  rate 

(radians/second) 

real 

right  ascension  (radians) 

real 

square  root  of  semi-major 

axis  (kilometers) 
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word  20  ADE 

word  21  IOCS 
word  22  TOES 
word  23  IWKNOS 
word  24  WS 

word  25  I  ED  ATS 
word  26  ADC 

If  ITYPE  equals  three  then  record 

word  1  TT 

word  2  ISAT(l) 

to  to 

wordNTR+1  ISAT(NTR) 

wordNTR+2      CR1(1) 

to  to 

wordNTR*2+l    CRl(NTR) 

wordNTR*2  +  2   CR2(1) 

to  to 

word  NTR*3+  1    CR2(XTR) 

wordNTR*3  +  2   DOPl(l) 

to  to 

wordNTR*4+l    DOPl(NTR) 

wordNTR*4  +  2   DOP2(l) 

to  to 

word  NTR*5  +  1    DOP2(NTR) 

word\TR*5  +  2     SGR1(1) 

to  to 

wordNTR*6+l     SGRl(NTR) 

wordXTR*6+2     SGR2(1) 


real 

real 
real 
intg 
real 


intg 
real 


age  of  data  ephemeris 

(seconds) 

clock  epoch  (seconds) 

ephemeris  epoch  (seconds) 

week  number 

argument  of  perigee 

(radians) 

sv  health 

age  of  data  clock  (seconds) 


two  will  contain  the  tracking  information. 

real  GPS  time-tag  for  data 

that  follows  (seconds) 
The  PRN  number  for  the 

intg  satellite  being  tracked 

on  trackers  i  to  NTR 


real 


real 


real 


real 


real 


The  LI  pseudorange  on 
each  tracker  (seconds) 

The  L2  pseudorange  on 
each  tracker 
(seconds) 

The  LI  continuous 
count  doppler  data  on 
each  tracker  (cycles) 

The  L2  continuous 
count  doppler  data  on 
each  tracker  (cycles) 

The  variance  on  the  LI 

pseudorange 

(kilometer*kilometcr) 

The  variance  on  the  L2 
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to  to  real 

wordNTR*7+l  SGR2(NTR) 

word  NTR*7  +  2  SGD  l(NTR) 

to  to  real 

wordNTR*8+l  SGDl(NTR) 

wordNTR*8  +  2  SGD2(1) 

to  to  real 

word  NTR*9  +  1  SGD2(NTR) 

wordNTR*9  +  2  MQVEL(l) 

to  to  intg 

word  NTR*10+  1  MQVEL(NTR) 


pseudorange 
(kilometer*kilometer) 

The  variance  on  the  LI 

doppler 

(kilometer*kilometer) 

The  variance  on  the  L2 

doppler 

(kilometer*kilometer) 

The  quality  vector  for 

all  trackers 

15  =  loss  of  lock 

14  =  bad  data,  no 

loss  of  lock 

1    =  good  range  data 

bad  doppler  data 

0    =  all  data  good 
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APPENDIX  D 
SOLTOM  SUBROUTINE  DESCRIPTIONS  AND  NESTING 


CARDIO  -  This  routine  reads  the  input  for  program  SOLTOM 

from  the  output  of  CON9TR. 

LLXYZ  -  This  routine  converts  position  in  the  geodetic 

coordinate  system  to  the  earth-fixed  geocentric 
coordinate  system. 

SIGXYZ  -  This  routine  converts  sigmas  in  the  geodetic 

system  into  sigmas  in  the  earth-fixed^ 
geocentric  coordinate  system. 

MATRIX  -  This  routine  performs  matrix  operations  by 

calling  routines  in  the  IMSL  library. 

LINV2F  -  This  routine  is  in  the  IMSL  and  does  a  matrix 

inversion. 

VMULFF  -  This  routine  is  in  the  IMSL  and  does  matrix 

multiplication. 

VML'LFM  -  This  routine  is  in  the  IMSL  and  does  a 

transpose  of  matrix  multiplication. 

INDAT  -  This  routine  will  pass  initial  values  to  some 

data  arrays. 

NEWSOL  -  This  routine  sets  up  the  configuration  for  a 

for  a  given  solution  routine. 

LSTSOL  -  This  routine  calculates  a  sequential  least 

squares  solution  for  the  receiver  position  and 
clock. 

GETDAT  -  This  routine  reads  the  next  set  of  satellite 

tracking  data  from  CON9TR  output  file  and  it 
will  also  update  weather,  Navigation  data, 
almanac  and  calibration. 

DEDPDT  -  This  routine  converts  the  accumulated  Doppler 

phase  data  into  range  difference  data  by 
correcting  phase  for"  offset  and  then 
differencing  it  with  data  from  before 
(if  before  data  is  good). 

DEFBR   -  This  routine  defines  bias  range 

CRDAT   -  This  routine  defines  correction  for  receiver 

bias. 

BSADR   -  This  routine  will  correct  data  for  receiver 

clock  bias,  drill  and  also  correct  GPS  time 
tag  for  the  end  of  week  cross  over. 

CRDAT2  -  This  routine  converts  data  into  range 

residuals  and  corrects  data  for  error  sources. 

ION      -  This  routine  corrects  for  ionospheric 

correction. 

IOXRD   -  This  routine  corrects  range  difference  data. 
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EPHDAT - 


CLOCK   - 

SATPOS  - 

ESTPOS  - 
ZANG     - 

RES      - 

SMOT    - 
TROP    - 
REL      - 
FSIG    - 

SOLTNS - 
SMPT1    - 

UDRRC   - 
SMPT2   - 

UDRDRC 
SMPT3   - 


This  routine  defines  the  following  values 
from  the  ephcmeris  data,  satellite  tracker 
number,  satellite  bias,  drift  and  aging  terms, 
semi-major  axis  of  satellite  orbit, 
eccentricity  of  the  satellite  orbit  and  epoch 
time  for  clock  terms. 

This  routine  adds  satellite  and  clock  error 
and  converts  it  to  range  measurements  in 
kilometers  and  also  performs  the  same  for  range 
difference  data. 

This  routine  will  find  the  position  of  the 
satellite  from  the  Navigation  Message  or  from 
the  Precise  Ephemeris. 

This  routine  finds  the  position  of  the 
satellite  from  the  Precise  Ephemeris  data. 

This  routine  finds  the  Zenith  angle  by  doting 
the  X,  Y,  Z  positions  of  the  satellite  receiver 

position  vector  to  the  vector  from  the  receiver 

o  the  satellite. 
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This  routine  finds  the  residual  by  calculating 
the  range  to  a  satellite  and  subtracting  it 
from  the  observed  ranse.   This  is  also 
performed  for  range  difference  data. 

This  routine  determines  an  earth  rotation 
correction. 

This  routine  finds  tropospheric  refraction 
correction. 

This  routine  finds  a  relativistic  correction 
to  the  data. 

This  routine  finds  the  variance  (sigma  squared) 
for  the  raw  range  and  range  difference  data  for 
the  satellite  being  tracked. 

This  routine  calls  the  different  routines  to 
sum  data  in  the  correct  part  of  the  "B"  matrix 
and  "E"  vector  depending  on  which  options  have 
been  set  by  the  user  for  type  of  solution. 

This  routine  calls  a  routine  to  define  additions 
to  "B"  matrix  and  "E"  vector  plus  directing  the 
called  routine  on  where  to  store  this 
information  for  the  "B"  matrix  and  "E"  vector 
(for  range  data  only). 

Update  ranee  solution  matrix  receiver  clocks. 
Then  store  data  in  matrix  so  that  solution  can 
computed. 

This  routine  calls  a  routine  to  define 
additions  to  "B"  matrix  and  "E"  vector  plus 
directing  the  called  routine  on  where  to  store 
this  information  for  "B"  matrix  and  "E"  vector 
(lor  range  difference  data  only) 

Update  range  difference  solution  matrix 
receiver  bias.    This  routine  stores  data  in 
matrix  so  a  solution  can  be  computed. 

This  routine  calls  a  routine  to  define 
additions  to  "B"  matrix  and  "E"  vector,  plus 
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NETRAP 

CROUT  - 
LINV3F  - 
POSAVG  ■ 

SMPT6   - 

UDRSC   - 

BELIM  - 
BMSTR  - 
SMPT7  - 

L'DRDSC 
SMPT8   - 

SMPT9   - 

UDCRSC 

ELIM     - 
SMPT10- 


dirccting  the  called  routine  on  where  to  store 
this  ini donation  for  the  "B"  matrix  and  "E" 
vector  (for  range  and  ranee  difference  data 
only). 

This  routine  uses  a  Newton-Raphson  method  to 
calculate  time  bias  adjustment  and  then  improve 
the  receiver  initial  position. 

This  routine  calls  an  IMSL  routine  that  finds 
the  determinate  of  a  matrix. 

This  is  an  IMSL  routine  that  calculates  the 
determinate  of  a  matrix. 

This  routine  finds  the  average  of  all  positions 
up  to  current  time  line. 

This  routine  calls  a  routine  to  define  additions 
to  "B"  matrix  and  "E"  vector  plus  directing  the 
called  routine  on  where  to  store  this 
information  (for  range  data  onlv,  satellite 
biases) 

This  routine  stores  data  in  matrix  so  solution 
so  solution  can  be  computed  after  the  ranee 
solution  matrix  and  satellite  clocks  have  been 
updated. 

This  routine  eliminates  certain  parameters  in 
the  "B"  matrix. 

This  routine  stores  data  in  a  matrix  so  a 
solution  can  be  computed. 

This  routine  calls  a  routine  to  define 
additions  to.  the  "B"  matrix  and  "E"  vector 
plus  directing  the  called  routine  on  where  to 
store  this  information  (for  range  difference 
data  only,  satellite  clock). 

This  routine  stores  data  in  a  matrix  so  that 
a  solution  can  be  computed  alter  updating  the 
range  solution  matrix  and  satellite  bias. 

This  routine  calls  a  routine  to  define  additions 
to  "B"  matrix  and  "E"  vector  plus  directing  the 
called  routine  where  to  store  this  information 
(for  range  and  range  difference  data  only). 

This  routine  calls  a  routine  to  define  additions 
to  the  "B"  matrix  and  "E"  vector  plus  directing 
the  called  routine  on  where  to  store  this 
information  (for  common  time  range  data  only). 

This  routine  stores  data  in  a  matrix  so  that 
a  solution  can  be  computed  after  updating  the 
common  range  solution  matrix  and  satellite 
clocks. 

This  routine  finds  data  needed  for  elimination 
and  then  calls  a  routine  to  do  that  elimination. 

This  routine  calls  a  routine  to  define 
additions  to  the  "B"  matrix  and  "E"  vector  plus 
directing  the  called  routine  on  where  to  store 
this  information  (for  Doppler  data  used  as  bias 
range  data). 
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BSRNG   -  This  routine  defines  additions  to  the  "B" 

matrix  and  "E"  vector. 

NEWITR  -  This  routine  resets  data  so  that  a  new  iteration 

can  be  done. 

BLOCK  DATA  -  This  routine  was  added  at  NPS  to  make  this 
program  compatible  with  the  FORI  RAN  77 
compilers  on  the  IBM-3033  mainframe 
computer  system.    It  contains  common 
statements  and  data  statements  that 
NPS's  compilers  would  not  permit  where 
they  were  placed  in  the  original  code  of 
program  SOLTOM. 


MAIN  PROGRAM  SOLTOM 
CARDIO 
LLHXYZ 
SIGXYZ 
MATRIX 

LINV2F 
VMULFF 
VMULFM 
IN  DAT 
NEWSOL 

LSTSOL 

MATRIX 
GETDAT 

DEDPDT 

DEFBR 
CRDAT 

BSADR 
CRDAT2 
ION 
IONRD 
LPIIDAT 
CLOCK 
SATPOS 

ESTPOS 
ZANG 
RES 
SMOT 
TROP 
REL 
FSIG 
SOLTNS 

SMPT1 

UDRRC 
SMPT2 

UDRDRC 
SMPT3 

UDRRC 
UDRDRC 
NET  RAP 

CROUT 

LINV3F 
CRDAT 
POSAVG 
SMPT6 

UDRSC 

BELIM 

MATRIX 

LINV2F 
VMULFF 
VMULFM 
BMSTR 
SMPT7 
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UDRDSC 

BELIM 
BMSTR 

SMPT8 

UDRSC 
UDRDSC 
SMPT9 

UDCRSC 

ELIM 

BELIM 
BMSTR 
SMPT10 

BSRNG 

CRDAT 
UDRSC 
NEWSOL 

NEWITR 

GETDAT 
BLOCK  DATA 
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APPENDIX  E 
INPUT  GUIDE  AND  SOLUTION  TYPES  FOR  PROGRAM  SOLTOM 


INPUT  GUIDE 

Var 

Type 

CARD  ONE 

col  1  to  10 
TSS 

real 

col  11  to  20 
TES 

real 

(  by  S.  L.  MeycrhofT,  NSWC) 


Description 

Format  (2 F  10.0,615) 
(This  line  can  oe  read  format  free) 

The  GPS  time  tag  for  data  at  start 
of  solution  (seconds) 


col  21  to  25  (right  justified) 
NPN  integer 

col  26  to  30  (right  justified) 
NPB  integer 


col  31  to  35  (right  justified) 
ITRY  integer 

col  36  to  40  (right  justified) 
NSAT  integer 


col  41  to  45  (right  justified) 
IGT  integer 


col  46  to  50  (right  justified) 
ISP  integer 


col  51  to  55  (right  justified) 
IAGE  integer 

col  56  to  60  (right  justified) 
ICLOCK  integer 


The  GPS  time-tag  of  last  data  point 
to  use  in  solution  (seconds) 


The  number  of  time  lines  of  data 
needed  before  first  solution 


The  number  of  time  lines  of  data 
between  solutions 


The  maximum  number  of  iterations 

The  maximum  number  of  trackers  to 
process  at  once  (range  1  to  4) 

=   1  if  data  source  is  Standard  9-track 

format 

=  2  if  data  source  GESAR  test  file 

format 

=  3  if  data  source  is  NAV  IF  tape 


=   1  for  satellite  position  from  NAV 

data 

=  2  for  satellite  position  from 

Celestial  Earth  fixed  trajectory 

files 


=   1  solve  for  clock  aging,  otherwise 
not  solved  for 


=  0  solve  clock  corrections  are  from 
NAV  data  message 
=   1  clock  corrections  are  post  fit 
data  from  cards 


CARD  TWO  {receiver  position) 
col  10 


Formal  {1 1 0,7 G  10.0) 
(  1  his  card  can  be  read  format  free) 
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I  FLAG 

integer 

col  11  to  20 
XLAT 

real 

col  21  to 
YLOG 

30 

real 

col  31  to 
ZHT 

40 

real 

col  41  to 
TB 

50 

real 

col  51  to 
TD 

60 

real 

col  61  to 
TO 

70 

real 

col  71  to 
AGE 

80 

real 

=  0  for  cartesian  coordinates 
=   1  for  geodetic  coordinates 

latitude  or  X  coordinate  for  position 
(degrees  or  kilometers) 

longitude  or  Y  coordinate  for  position 
(degrees  or  kilometers) 

height  or  Z  coordinate  for  position 
(kilometers) 

time  bias  for  receiver  clock 
(kilometers) 

time  drift  for  receiver  clock 
(kilometers/second) 

epoch  for  clock  drift  (given  and 
solution)  If  negative,  then  first 
used  time  line  is  the  epoch  for  the 
solution  clock  drift 


aging  for  receiver  clock 
(kilometers/ second*  second) 


CARD  THREE  {Reference  Geoid)  Format  (2 G  10.0) 

(This  card  can  be  read  format  free) 


semi-major  axis  of  the  reference  geoid 
(kilometers) 

Obliquity  of  the  reference  geoid 

Position)  Format  (7G  10.0) 

(This  card  can  be  read  format  free) 

sigma  of  X  or  latitude  coordinate 
(kilometers) 

sigma  of  Y  or  longitude  coordinate 
(kilometers) 

sigma  of  Z  or  height  coordinate 
(kilometers) 

sigma  of  refraction  scaling  constant 

sigma  on  receiver  or  satellite  1  time 
bias  (kilometers) 


sigma  on  receiver  or  satellite  1 
clock  drift  (kilometers/second) 


col  1  to  10 
AG 

real 

col  11  to  20 
ROBL 

real 

CARD  FOUR  (Si 

'gmas 

col  1  to  10 
S1GLT 

real 

col  1 1  to  20 
SIGLG 

real 

col  21  to  30 
SIGHT 

real 

col  31  to  40 
SIG(l) 

real 

col  41  to  50 
SIG(2) 

real 

col  51  to  60 
SIG(3) 

real 
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col  61  to  70 
SIG(4) 


real 


sigma  for  aging  term  or  satellite  2 
bias  (kilometers/second*sccond) 


CARD  FIVE  ( Sigmas  for  Position)  Format  (7G 1 0. 0) 

(This  card  can  be  read  format  free) 


col  1  to  10 
SIG(5) 

real 

siema  on  satellite  2  drift 
(kilometers/ second) 

col  11  to  20 
SIG(6) 

real 

sigma  for  satellite  3  bias 
(kilometers) 

col  21  to  30 
SIG(7) 

real 

sigma  for  satellite  3  drift 
(kilometers/second) 

col  31  to  40 
SIG(8) 

real 

sigma  for  satellite  4  bias 
(kilometers) 

col  41  to  50 
SIG(9) 

real 

sigma  for  satellite  4  drift 
(kilometers) 

SOLUTION  TYPES 

CARD  SIX  (Types  of  Solutions  to  do)  Format  (1615) 

(This  card  cannot  be  read  format  free) 


=   1  for  range  solution  receiver  clock 
0  =  no 


=   1  for  range  difference  solution 
receiver  clock 
=  0  no 


=   1  for  range  and  range  difference 
solution  receiver  clock 
=  Ono 


=   1  for  common  time  range  solution 
=  0  no 


=   1  for  Newton- Raphson  range  solution 
=  0  no 


=   1  for  range  solution  satellite 
=  0  no 

=   1  for  range  difference  solution 
=  0  no 


col  5 
NSOL(l) 

integer 

col  10 
NSOL(2) 

integer 

col  15 
NSOL(3) 

integer 

col  20 
NSOL(4) 

integer 

receiver  clock 

col  25 
NSOL(5) 

integer 

col  30 
NSOL(6) 

integer 

clocks 

col  35 
NSOL(7) 

integer 

satellite  clock 

col  40 
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NS0L(8) 


col  45 
NSOL(9) 


col  50 
XSOL(IO) 


integer 


1  for  range  and  range  difference 


solution  satellite  clock 

=  0  no 


integer 
satellite  clock 

integer 
range  data 


=  1  for  common  time  range  solution, 

=  0  no 

=  1  for  doppler  data  used  as  bias 

=  Ono 


CARD  SEVEN  ( Variance  Bias  Factor)  Format  (2G10.0) 

(This  card  can  be  read  format  free) 

col  1  to  10 
K2 


col  11  to  20 

K4 


real 


real 


bias  factor  for  variance  on  range 
data 


bias  factor  for  variance  on  doppler 
data 


CARD  EIGHT  (Print  Options)  Format  (515) 

(7  his  card  can  be  read  format  free) 


col  5 
IPRT(l) 

col  10 
IPRT(2) 


col  15 
IPRT(3) 

col  20 
IPRT(4) 


col  25 
IPRT(5) 

col  30 
IPRT(6) 


col  35 
IPRT(7) 


col  40 
IPRT(8) 


integer 

integer 
terms 

integer 

integer 
data  sigma 

integer 

integer 
corrections 

integer 

lat,  long,  height 

integer 


=   1  print  position  solution  (X,Y,Z) 
=  0  do  not  print 

=   1  print  solution  sigmas  and  clock 

=  0  do  not  print 

=   1  print  data  and  time  tags 
=  0  do  not  print 

=   1  print  residuals,  range,  trop, 

=  0  do  not  print 

=   1  print  satellite  position 
=  0  do  not  print 

=   1  print  residuals  and  all  data 

=  0  do  not  print 

=   1  print  solution  position  in 

=  0  do  not  print 

=   1  print  real  time  solutions  from 
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9-track 

=  0  do  not  print 

CARD  NINE  (Weather  Data.  Elev.  Cutoff)  Formal  (4GW.0) 

( This  card  can  be  read  format  free) 

col  1  to  10 
PB  real  pressure  in  millibars 

col  1 1  to  20 
TC  real  temperature  in  degrees  Centigrade 

col  21  to  30 
HD  real  relative  humidity  in  percent 

col  31  to  40 
EL  real  elevation  angle  cutolT  (degrees) 

CARD  TEN  (Station  Number)  Format  (15) 

(This  card  can  be  read  format  free) 

col  1  to  5 
ISTA  integer  station  identification 

CARD  ELEVEN  (Tracker  1  Clock  Parameters)  Format  (3G  10.0) 

col  1  to  10 
AO  real  clock  bias  of  satellite  tracked  on 

tracker  1 

col  1 1  to  20 
Al  real  clock  drift  of  satellite  tracked  on 

tracker  1 

col  21  to  30 
A2  real  clock  aging  of  satellite  tracked  on 

tracker! 

CARD  TWELVE  (Tracker  2  Clock  Parameters)  Format  (3G10.0) 

(This  card  can  be  read  lormat  free) 

col  1  to  10 
AO  real  clock  bias  of  satellite  tracked  on 

tracker  2 

col  1 1  to  20 
Al  real  clock  drift  of  satellite  tracked  on 

tracker  2 

col  21  to  30 
A2  real  clock  aging  of  satellite  tracked  on 

tracker  7 

CARD  THIRTEEN  (Tracker  3  Clock  Parameters)  Format  (3G10.0) 

(This  card  can  be  read  format  free) 

col  1  to  10 
AO  real  clock  bias  of  satellite  tracked  on 

tracker  3 

col  1 1  to  20 
Al  real  clock  drift  of  satellite  tracked  on 

tracker  3 

col  21  to  30 
A2  real  clock  aging  of  satellite  tracked  on 

tracker  J 
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CARD  FOURTEEN  {Tracker  4  Clock  Parameters)  Format{3G  10.0) 

( I  his  card  can  be  read  format  irce) 

clock  bias  of  satellite  tracked  on 
tracker  4 

clock  drift  of  satellite  tracked  on 
tracker  4 

clock  aging  of  satellite  tracked  on 
tracker  4 


col  1  to  10 
AO 

real 

col  11  to  20 
Al 

real 

col  21  to  30 
A2 

real 
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APPENDIX  F 
KALMN2  SUBROUTINE  DESCRIPTIONS  AND  NESTING 


CARDIO  -    .  This  routine  reads  input  for  program. 

LLHXYZ  -  This  routine  converts  latitude,  longitude  and 

height  to  X,  Y.Z  in  earth-fixed  geocentric 
cartesian  coordinates. 

SIGXYZ  -  This  routine  converts  sigmas  in  the  geodetic 

coordinate  system  into  sigmas  in  the  earth- 
fixed  geocentric  coordinate  system. 

MATRIX  -  This  routine  performs  matrix  operations  bv 

calling  the  proper  routine  in  the  IMSL  library. 

LINV2F  -  This  routine  is  in  the  IMSL  and  performs  a 

matrix  inversion. 

VMULFF  -  This  routine  is  in  the  IMSL  and  performs  a 

matrix  multiplication. 

VMULFM  -  This  routine  performs  a  transpose  of 

multiplication  of  matrices. 

INDAT  -  This  routine  gives  the  initial  values  to  the 

"B"  matrix  (L  in  Kalman  filter) 

TIMEUP  -  This  routine  uses  a  UDU'  factorization  method 

to  update  state  matrix. 

PROMPT  -  This  routine  aids  in  the  UDU'  factorization. 

GETDAT  -  This  routine  reads  the  next  set  of  satellite 

tracking  data  from  CON9TR  output  file. 

NEWRAN  -  This  routine  reduces  noise  and  multipath  in  the 

range  data. 

DEDPDT  -  This  routine  converts  Doppler  phase  into  range 

differenced  data. 

DEFBR  -  This  routine  defines  bias  range. 

CRD  AT   -  This  routine  defines  corrections  for  receiver 

bias. 

BSADR   -  This  routine  will  correct  data  for  receiver 

clock  bias  and  drift  and  also  correct  GPS  time- 
tag  for  the  end  of  week  cross  over. 

CRDAT2  -  This  routine  converts  data  into  range  residuals 

and  corrects  data  for  error  sources. 

ION      -  This  routine  corrects  for  ionospheric 

refraction. 

IONRD   -  This  routine  corrects  range  dilfcrcnce  data 

EPIIDAT  -  This  routine  defines  the  values  from  the 

cphemeris  data. 

CLOCK    -  This  routine  adds  satellite  and  clock  error  and 
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SATPOS - 

ESTPOS  - 
ZANG     - 

RES      - 
SMOT    - 

TROP    - 

REL      - 

FSIG     - 

GDOP    - 
LINV3F  - 

KALFIL- 

INKAL   - 

UDFACT 

RKAL    - 

UPDATE ■ 

DPKAL   - 

COXSTR  - 

XYZLLII  • 

PRNSOL  - 
COVEVL - 


converts  it  to  range  measurements  in  kilometers 
and  also  perform  the  same  for  range  difference 
data. 

This  routine  will  find  the  position  of  the 
satellite  from  the  navigation  data  message  or 
from  the  Celestial  Earth- Fixed  Trajectory 
(Precise  Ephemeris). 

This  routine  finds  the  position  of  the 
satellite  from  the  Precise  Ephemeris  data. 

This  routine  finds  the  zenith  angle  by  doting 
the  X,  Y,  Z  positions  of  the  satellite  "receiver 
position  vector  to  the  vector  from  the  receiver 
to  the  satellite. 

This  routine  finds  the  residual  for  ranges. 

This  routine  determines  an  earth  rotation 
correction. 

This  routine  calculates  a  tropospheric 
refraction  correction. 

This  routine  finds  a  relativistic  correction 
to  the  data. 

This  routine  calculates  the  variance  for  raw 
pseudorange  and  range  differenced  data  for  the 
satellite  being  tracked. 

This  routine  if  turned  on,  writes  a  GDOP  file. 

This  routine  is  in  the  IMSL  and  performs  a 
matrix  inversion 

This  routine  uses  corrected  data  to  do  the 
Kalman  filter. 

This  routine  sets  up  initial  covariance  matrix 
in  UDU'  format. 

This  routine  puts  Apriori  covariance  matrix 
into  UDU'  factorized  matrix. 

This  routine  calculates  partial  derivatives  of 
the  data  equation  for  pseudorange  data. 

This  routine  will  update  the  Kalman  solution 
for  one  data  point. 

This  routine  calculates  partial  derivatives  of 
the  data  equation  for  delta  (Doppler)  range 
differences. 

This  routine  will  constrain. the  solution  to 
the  surface  of  a  sphere  (height)  and  constrain 
velocity. 

Converts  X,  Y,  Z  in  an  earth-fixed  geocentric 
coordinate  system  to  geodetic  coordinate  system. 

This  routine  prints  solutions. 

This  routine  finds  the  covariance  matrix  from 
the  UDU"  factorized  matrix  in  U. 


BLOCK  DATA  -  This  routine  was  added  at  NPS  to  handle 


67 


data  statements  and  common  statements  that 

had  to  be  moved  in  the  code  or  written,  to 

be  compatible  with  the  1BM-3033  FORTRAN  77 

compilers. 


MAIN  PROGRAM  KALMN2 
CARDIO 

LLHXYZ 
SIGXYZ 

MATRIX 

LINV2F 
VMULFF 
VMULFM 
INDAT 

TIMEUP 

PROMPT 
GETDAT 

DEDPDT 

DEFBR 
CRDAT 

BSADR 
CRDAT2 

NEW  RAN 

IONRD 
ION 
IONRD 
EPHDAT 
CLOCK 2 
CLOCK 
SATPOS 

ESTPOS 
ZANG 
RES 
SMOT 
TROP 
REL 
FSIG 
GDOP 

LINV3F 
K  A  LITE 

INKAL 

UDFACT 
LLHXYZ 
RGKAL 

UPDATE 
DPKAL 

UPDATE 
CONST  R 

UPDATE 
XYZLLII 
PRNSOL 

COVEVL 
BLOCK  DATA 
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Var         Type 

CARD  ONE 

col  1  to  10 
TSS 


col  1 1  to  20 
TES 


APPENDIX  G 
INPUT  GUIDE  FOR  PROGRAM  KALMN2 

(  by  S.  L.  MeyerhofT,  NSWC) 

Description 

Format  (2F10.0.6I5) 
(This  card  can  be  read  format  free) 


real 


real 


col  21  to  25  (right  justified) 
NPN  integer 


col  26  to  30  (right  justified) 
NPB  integer 


col  31  to  35  (rieht  justified) 
ITRY  integer 


col  36  to  40  (right  justified) 
NSAT  integer 


col  41  to  45  (right  justified) 
IGT  integer 


col  46  to  50  (right  justified) 
ISP  integer 


col  51  to  55  (right  justified) 
ICLOCK  integer 


col  56  to  60  (right  justified) 
1AGE  integer 


col  61  to  65  (right  justified) 
IAGE  integer 


The  GPS  time-tag  for  data  at  start 
of  solution  (seconds) 

The  GPS  time-tag  oflast  data  point 
to  use  in  solution  (seconds) 


The  number  of  time  lines  between 
solutions 


The  number  of  solutions  between 
solutions  being  printed 

The  number  of  solutions  printed 
covariance  matrix  prints 


The  number  of  trackers  to  process 
in  the  solution  (range  1  to  4) 


=   1  if  data  source  is  standard 

9-track  format 

=  2  if  data  source  is  other 


=   1  for  satellite  position  from 
NAVIGATION  DATA  MESSAGE 

=  2  for  satellite  position  from 
Celestial  Earth  Fixed 
Trajectory  files 


=  0  clock  corrections  are  from 
NAVIGATION  DATA  MESSAGE 
=   1  clock  corrections  are  post 
fit  data  from  cards 


=   1  make  GDOP  file 
otherwise  no  file  is  made 


=   1  solve  for  clock  aging 
otherwise  no  file  is  made 


CARD  TWO  (Receiver  Position)      Format  {IIOJGJO.O) 

(This  card  can  be  read  format  free) 

col  10 
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I  FLAG 


integer  =  0  for  Cartesian  coordinates 

=   1  for  Geodetic  coordinates 


col  11  to  20 
XL  AT 


col  21  to  30 
YLOG 


col  31  to  40 
ZHT 


col  41  to  50 
IB 


col  51  to  60 
YD 


col  61  to  70 
TO 


col  71  to  80 
AGE 


real  latitude  or  X  coordinates  for 

position  (degrees  or  kilometers) 

real  longitude  or  Y  coordinate  for 

position  (degrees  or  kilometers) 

real  height  or  Z  coordinate  for 

position  (kilometers) 

real  time  bias  for  receiver  clock 

(kilometers) 

real  time  drift  for  receiver  clock 

(kilometers/second) 

real  epoch  for. clock  drift  (given  and 

solution)  if  negative,  then  first 
used  time  line  is  the  epoch  for 
the  solution  clock  drift 


real  aging  for  receiver  clock 

(kilometers/  second*second) 

CARD  THREE  {Reference  Geoid)        Format  (2G10.0) 

(This  card  can  be  read  format  free) 
col  1  to  10 
AG  real  semi-major  axis  of  the  reference 

geoid  (kilometers) 


obliquity  of  the  reference  geoid 

Position)     Format  (7 G  10.0) 

(This  card  can  be  read  format  free) 

sigma  of  X  or  latitude  coordinate 
(kilometer) 

sigma  of  Y  or  longitude  coordinate 
(kilometer) 

sigma  of  Z  or  height  coordinate 
(kilometers) 

sigma  on  refraction  scaling  constant 


sigma  on  receiver  time  bias 
(kilometer) 


sigma  on  receiver  clock  drift 
(kilometer/ second) 


col  11  to  20 
ROBL 

real 

CARD  FOUR  (Sigmasfc 

col  1  to  10 
SIGLT 

real 

col  1 1  to  20 
SIGLG 

real 

col  21  to  30 
SIGHT 

real 

col  31  to  40 
SIG(l) 

real 

col  41  to  50 
SIG(2) 

real 

col  51  to  60 
SIG(3) 

real 

col  61  to  70 

70 


SIG(4) 


real 


sigma  of  velocity  in  X  direction 
(kilometers/ second) 


CARD  FIVE  (Sigmasfor  Velocity)      format  (7GI0.0) 

(This  card  can  be  read  format  free) 


col  1  to  10 
SIG(5) 

col  1 1  to  20 
SIG(6) 

col  21  to  30 

SIG(7) 

col  31  to  40 
SIG(8) 

col  41  to  50 
SIG(9) 


real 
real 
real 
real 
real 


sigma  of  velocity  in  Y  direction 
(kilometers/second) 

sigma  of  velocity  in  Z  direction 
(kilometers/ second) 

sigma  on  receiver  time  bias 
(kilometers) 

sigma  on  receiver  clock  drift 
(kilometers/second) 

sigma  not  used 


CARD  SIX  (Types  of  solutions  to  do)     Format  (1615) 

(This  card  cannot  be  read  format  free) 


col  5 
NSOL(l) 

col  10 
NSOL(2) 


col  15 
NSOL(3) 

col  20 
NSOL(4) 


col  25 
NSOL(5) 


col  30 
NSOL(6) 


col  35 
XSOL(7) 


col  40 
NSOL(8) 


integer  =   1  for  range  position  solution 

=  0  no 


integer  =   1  Height  constraint  on  position 

=  0  constraint  only  if  less  than  4 
satellites 


integer  =   1  write  position  solution  on  file 

=  0  no  solution  saved 


integer  =   1  to.  solve  for  refraction 

correction  in  position 
solution 
=  0  no 


integer  =   1  to  solve  for  refraction 

correction  in  velocity 
solution 
=  0  no 


integer  =   1  to  solve  for  clock  drift 

from  position  solution 
=  0  no 


integer  =   1  do  a  velocity  solution 

(this  has  not  been  used 
before) 
=  0  no 


integer  =   1  do  a  time  update  of 

covariance  matrix 
=  0  no 
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col  45 
NSOL(9) 

col  50 
NSOL(IO) 


integer 
integer 


=   1  not  used 
=  Ono 


=   1  for  a  constraint  on 
velocity  solution 

=  Ono 


CARD  SEVEN  (Variance  Bias  Factor)   Format  (2G10.0) 

(This  card  can  be  read  format  free) 


col  1  to  10 
K2 


col  1 1  to  20 
K4 


real 


real 


bias  factor  for  variance  on  range 
data 


bias  factor  for  variance  on 
Doppler  data 


CARD  EIGHT  (Print  Options)      Format  (515) 

(This  card  cannot  be  read  format  free) 

=   1  print  position  solution 

(X.YS) 

=  0  do  not  print 

=   1  print  solution  sigmas  and 

clock  terms 

=  0  do  not  print 

=   1  print  data  and  time  tags 
=  0  do  not  print 

=   1  residuals,  range,  trop, 

data  sigma 

=  0  do  not  print 

=   1  print  satellite  position 
=  0  do  not  print 

=   1  residuals  and  all  data 

corrections 

=  0  do  not  print 

=   1  print  position  in  lat,  long, 
and  height 

=   1  print  real  time  solutions  from 
9-track 

CARD  NINE  (weather  Data,  Elev.  Cutoff)  Format  (46 10.0) 

( 1  his  card  can  be  read  lormat  free) 


col  5 
IPRT(l) 

integer 

col  10 
IPRT(2) 

integer 

col  15 
IPRT(3) 

integer 

col  20 
IPRT(4) 

integer 

col  25 
IPRT(5) 

integer 

col  30 
IPRT(6) 

integer 

col  35 
1PRT(7) 

integer 

col  40 
IPRT(8) 

integer 

col    1  to  10 
PB 

col  1 1  to  20 
TC 


real 
real 


pressure  in  millibars 

temperature  in  degrees  Centigrade 
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col  21  to  30 
HD 

real 

col  31  to  40 
EL 

real 

CARD  TEN  (Station  Nu 

col    1  to  10 

Q(l) 

real 

col  11  to  20 
0(2) 

real 

col  21  to  30 
Q(3) 

real 

col  31  to  40 
0(4) 

real 

col  41  to  50 
SIGPS 

real 

col  51  to  60 
SIGVL 

real 

CARD  ELEVEN  (Diago. 

col    1  to  10 
PHI(1,1) 

real 

col  11  to  20 
PHI(2,2) 

real 

col  21  to  30 
PHI(3,3) 

real 

col  31  to  40 
PH  1(4,4) 

real 

col  41  to  50 
PI  11(5,5) 

real 

col  51  to  60 
PHI(6,6) 

real 

col  61  to  70 
PHI(7,7) 

real 

relative  humidity  in  percent 

elevation  angle  cutofT  (degrees) 

xber)  E or  mat  (Free) 

(This  card  can  be  read  format  free) 

time  propagation  noise  for  X 
velocity  (kilometers) 

time  propagation  noise  for  Y 
velocity  (kilometers) 

time  propagation  noise  for  Z 
velocity  (kilometers) 


time  propagation  noise  for  time 
drift 


sigma  on  position  constraint 


sigma  on  velocity  constraint 

xl  elements  of  state  transition  matrix) 
(This  card  can  be  read  format  free) 

state  transition  matrix  for 
parameter  1 


state  transition  matrix  for 
parameter  2 

state  transition  matrix  for 
parameter  3 

state  transition  matrix  for 
parameter  4 

state  transition  matrix  for 
parameter  5 

state  transition  matrix  for 
parameter  6 

state  transition  matrix  for 
parameter  7 

CARD  TWELVE  (Diagonal  elements  of  stale  transition  matrix) 

(This  card  can  be  read  lormat  Ircc) 


col    1  to  10 
PHI(8,8) 

col  1 1  to  20 


real 


state  transition  matrix  for 
parameter  8 
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PH  1(9,9)  real  state  transition  matrix  for 

parameter  9 

CARD  THIRTEEN  (Height  for  position  constraint) 

('I  "his  card  can  be  read  format  free) 

col    1  to  10 
II T  real  Height  above  reference 

ellipsoid  to  constrain 
position  solution  to 
(kilometer) 
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APPENDIX  H 
ABBREVIATIONS  AND  ACRONYMS 


GPS 

NTS 

DMA 

NOAA 

NSWC 

sv 

SEP 

P-code 

C/A-code 

WGS 

PDOP 

HDOP 

VDOP 

TDOP 

SEF 

CIO 

ARL 

IBM 
CDC 

m 

km 

s 

Hz 

MHz 


Global  Positioning  System 

Naval  Postgraduate  School  Monterey  CA 

Defense  Mapping  Agency 

National  Oceanic  and  Atmospheric  Administration 

Naval  Surface  Weapons  Center 

Space  Vehicle 

Spherical  Error  Probability 

Precision  code 

Coarse  Acquisition  code 

World  Geodetic  System 

Position  Dilution  of  Precision 

Horizontal  Dilution  of  Precision 

Vertical  Dilution  of  Precision 

Time  Dilution  of  Precision 

Standardized  Exchange  Format 

Conventional  International  Origin 

Applied  Research  Laboratories  at  the  University  of  Texas 
at  Austin 

International  Business  Machines 

Computer  Data  Corporation 

meters 

kilometers 

second 

Hertz 

Megahertz 
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